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Executive  Summary 

This  report  summarizes  progress  for  the  first  two  quarters  of  SERDP  project  MR-2226. 
The  primary  objective  of  the  initial  phase  of  the  project  is  the  development  of  a  decision 
support  system  (DSS)  to  predict  classification  performance  given  environmental  variables. 
This  report  details  progress  on  the  following  tasks: 

(f)  Prediction  of  target  detection  thresholds.  Current  methods  for  setting  target  detection 
thresholds  are  based  on  the  worst  case  (i.e.  minimum)  amplitude  response  for  a  given  target 
of  interest.  The  geophysical  system  verification  (GSV)  calculator  previously  developed  by 
ESTCP  assumes  that  this  worst  case  corresponds  to  a  horizontal  target  oriented  cross-track. 
We  show  that,  depending  on  sensor  geometry  and  offset  from  the  center  of  the  array,  the 
minimum  target  response  for  a  horizontal  target  is  not  always  in  a  cross  track  azimuthal 
orientation.  In  addition,  we  show  that  the  worst  case  target  amplitude  can,  in  practice, 
produce  a  prohibitively  large  number  of  target  picks.  To  reduce  the  number  of  picks,  we 
simulate  the  distribution  of  target  response  amplitudes  assuming  a  uniform  distribution  of 
target  azimuth  and  location  across  the  sensor  swath.  We  then  show  that  setting  a  threshold 
that  provides  a  very  high  confidence  of  detecting  targets  of  interest  (e.g.  99  %)  can  drastically 
reduce  the  number  of  target  picks. 

(2)  Noise  estimation.  Most  synthetic  analyses  use  discrete  background  measurements  to 
characterize  noise  at  a  site.  However,  we  find  that  a  simple  noise  model  derived  solely  from 
background  measurements  (e.g.  independent,  Gaussian  noise)  cannot  accurately  reproduce 
the  distributions  of  dipole  polarizabilities  recovered  from  field  data.  We  instead  use  the 
recovered  polarizability  distribution  of  industry  standard  objects  (ISOs)  as  a  barometer 
for  characterizing  noise  across  a  site.  We  develop  techniques  to  invert  the  covariance  of 
polarizabilities  within  the  ISO  class  to  recover  an  effective  covariance  of  the  noise  on  the  data 
across  the  site.  The  resulting  noise  is  highly  correlated  and  only  approximately  Gaussian. 

(3)  Performance  prediction.  Our  aim  here  is  to  predict  the  receiver  operating  character¬ 
istic  (ROC)  that  would  be  obtained  for  a  sensor  and  targets  under  specified  environmental 
conditions.  For  a  given  site,  we  associate  the  noise  covariance  derived  from  ISO  polariz¬ 
abilities  with  the  median  location  of  ISO  items,  relative  to  the  sensor.  With  this  noise 
covariance,  we  can  predict  the  covariance  of  target  polarizabilities  at  any  other  location  via 
a  straightforward  linear  analysis. 

Next,  we  derive  analytic  expressions  for  the  misfit  of  log-transformed  polarizabilities  with 
respect  to  reference,  or  library,  polarizabilities.  This  misfit  is  the  statistic  most  commonly 
used  to  discriminate  between  targets  of  interest  (TOI)  and  non-TOI  using  cued  sensor  data. 
The  distribution  of  the  decision  statistic  then  maps  to  the  predicted  ROC.  This  machinery 
allows  us  to  efficiently  predict  how  the  spatial  distribution  of  targets  affects  classification 
performance.  We  derive  analytic  expressions  for  the  covariance  of  target  polarizabilities 
under  a  uniform  distribution  of  target  orientations.  We  show  that  target  dip  has  an  impor¬ 
tant  effect  on  classification  performance,  with  a  uniform  distribution  of  dips  producing  ROC 
performance  that  is  intermediate  to  horizontal  and  vertical  cases. 

Finally,  we  extend  our  performance  prediction  methods  to  dynamic  sensors  and  develop 
techniques  to  account  for  the  effect  of  varying  data  density  on  classification  performance. 
We  derive  expressions  for  the  distribution  of  the  polarizability  decay  parameter,  which  is 
typically  used  to  classify  data  acquired  with  the  EM-61  sensor.  This  provides  a  means  to 
predict  the  baseline  ROC  for  a  site,  and  so  to  predict  whether  cued  interrogation  will  produce 
a  significant  improvement  in  classification  performance. 


1 


Contents 

Executive  Summary  2 

List  of  Figures  2 

List  of  Tables  3 

List  of  Acronyms  4 

1.  Introduction  5 

2.  Background  5 

2.1.  Data  acquisition  5 

2.2.  Feature  extraction  6 

2.3.  Classification  9 

3.  Predicting  target  response  thresholds  11 

4.  Performance  prediction  16 

4.1.  Performance  prediction  for  cued  interrogation  17 

4.1.1.  Noise  estimation  17 

4.1.2.  Predicting  the  distributions  of  polarizabilities  27 

4.1.3.  Predicting  the  ROC  35 

4.2.  Performance  prediction  for  dynamic  sensors  42 

5.  Conclusions  and  further  work  48 

References  49 

Appendix  A.  Summary  of  ESTCP  demonstration  data  sets  51 

Appendix  B.  Expected  polarizability  covariance  under  a  uniform  distribution  of 

target  orientations  53 


2 


List  of  Figures 

1  Electromagnetic  induction  (EMI)  survey  5 

2  Sequential  inversion  approach  9 

3  Computing  the  ROC  from  the  score  distributions  of  TOI  and  non-TOI  11 

4  Anomaly  amplitude  analysis  for  a  three  coil  EM-61  array  12 

5  Anomaly  amplitude  analysis  for  an  asymmetric  three  coil  EM-61  array.  13 

6  Cumulative  distribution  of  response  amplitudes  14 

7  Cumulative  distribution  of  response  amplitudes  at  JDR  15 

8  Distribution  of  response  amplitudes  for  targets  distributed  uniformly  in  azimuth 

and  dip.  15 

9  Methods  for  polarizability  prediction  17 

10  Background  noise  standard  deviations  for  MetalMapper  sensor  at  Camp  Butner  and 

Pole  Mountain.  18 

11  Components  of  the  observed  secondary  field  for  MetalMapper  dynamic  data  acquired 

at  Camp  Butner  19 

12  Estimated  ISO  polarizabilities  at  Camp  Beale  and  Pole  Mountain.  19 

13  Synthetic  example  of  estimation  of  effective  noise  standard  deviation  21 

14  Estimates  of  effective  noise  standard  deviations  for  MetalMapper  Beale  Parsons  and 

Pole  Mountain  data  sets  22 

15  Comparison  of  actual  and  synthetic  estimated  ISO  polarizabilities  at  Pole  Mountain  23 

16  Example  fit  of  a  sum  of  exponential  decays  to  estimated  polarizabilities  for  an  ISO 

target  24 

17  MetalMapper  data  predicted  using  unsmoothed  and  smoothed  polarizability 

estimates  25 

18  Fits  to  noise  distribution  derived  from  ISO  target  polarizabilities,  for  channel  42  of 

MetalMapper  data  at  Pole  Mountain  25 

19  Estimated  normal  mixture  model  parameters  at  each  MetalMapper  channel  26 

20  Structure  of  the  noise  covariance  matrices  at  Camp  Beale  and  Pole  Mountain  28 

21  Effect  of  target  depth  on  polarizability  uncertainty  29 

22  Dependence  of  polarizability  standard  deviations  on  horizontal  target  location  30 

23  Predicted  ISO  polarizabilities  for  four  spatial  distribution  scenarios  31 

24  ISO  polarizabilities  for  TEMTADS2x2  and  BUDhh  handheld  sensors  at  Camp  Beale  32 

25  Standard  deviation  of  positional  error  as  a  function  of  target  position  34 

26  Targets  used  for  ROC  prediction  37 

27  ROC  prediction  for  the  MetalMapper,  Pole  Mountain  noise  37 

28  ROC  prediction  for  the  MetalMapper,  Camp  Beale  noise.  38 

29  Effect  of  target  orientation  on  decision  statistic  and  ROC  39 

30  Numerical  simulations  of  primary  polarizability  estimation  for  vertical  and  horizontal 

ISO  targets.  40 

31  Total  polarizabilities  estimated  for  ISO  targets  at  Camp  Beale.  42 

32  EM-61  size-decay  feature  space,  Camp  Beale  handheld  area.  43 


3 


33  Dependence  of  EM61  total  polarizability  standard  deviations  at  Camp  Beale  on 

cross  track  spacing  44 

34  Dependence  of  lognormal  decay  distribution  on  target  size  and  correlation  coefficient  45 

35  Performance  prediction  for  the  EM61,  Camp  Beale  noise.  47 

List  of  Tables 

1  EM  sensors  used  for  UXO  classification  7 


4 


List  of  Acronyms 


BUD 

Berkeley  UXO  Discriminator 

BUDHH 

Handheld  Berkeley  UXO  Discriminator 

CDF 

Cumulative  Density  Function 

cm 

centimeter 

DSS 

Decision  Support  System 

EM 

Electromagnetic 

EMI 

Electromagnetic  induction 

ESTCP 

Environmental  Security  Technology  Certification  Program 

FAR 

False  Alarm  Rate 

FPF 

False  Positive  Fraction 

FLBGR 

Former  Lowry  Bombing  and  Gunnery  Range 

GSV 

Geophysical  System  Verification 

ISO 

Industry  Standard  Object 

JDR 

Jeep  and  Demolition  Range 

m 

Meter 

mm 

Millimeter 

MPV 

Man-Portable  Vector  Sensor 

MSEMS 

Man-Portable  Simultaneous  EMI  and  Magnetometer  System 

ms 

Millisececond 

MTADS 

Multi-Sensor  Towed  Array  Detection  System 

mV 

MilliVolt 

PI 

Principal  Investigator 

POC 

Point  of  Contact 

QC 

Quality  Control 

ROC 

Receiver  Operating  Characteristic 

SLO 

San  Luis  Obispo 

SNR 

Signal  to  Noise  Ratio 

TEM 

Time-domain  electromagnetic 

TEMTADS 

Time  Domain  Electromagnetic  Towed  Array  Detection  System 

TOI 

Target  of  Interest 

TPF 

True  Positive  Fraction 

UBC-GIF 

University  of  British  Columbia  Geophysical  Inversion  Facility 

UXO 

Unexploded  Ordnance 

5 


1.  Introduction 

While  significant  advances  have  been  made  in  the  acquisition  and  processing  of  geophysical 
data  for  classification  of  buried  munitions,  the  success  of  any  classification  strategy  strongly 
depends  upon  site  characteristics,  including  range  of  munitions  types  and  clutter,  geological 
background,  topography  and  vegetation.  In  this  project  our  goal  is  to  develop  and  validate 
the  components  of  a  decision  support  system  (DSS)  that  will  help  site-managers  and  teams 
design  surveys  and  data  processing  strategies  to  achieve  optimal  classification  performance 
at  the  lowest  attainable  cost  for  a  given  site. 

2.  Background 

Advanced  classification  of  buried  munitions  requires  a  number  of  steps: 

(1)  Data  acquisition:  detection  of  buried  targets  with  a  geophysical  sensor. 

(2)  Feature  extraction:  characterization  of  each  target  with  features  estimated  through 
inversion  of  a  parameterized  physics  based  forward  model. 

(3)  Classification:  prioritization  of  targets  for  digging  using  estimated  features. 

The  ultimate  goal  of  this  processing  is  to  identify  all  targets  of  interest  with  a  minimal 
number  of  false  alarms. 

2.1.  Data  acquisition.  In  the  data  acquisition  stage,  a  geophysical  sensor  is  deployed 
at  the  site.  Time-domain  electromagnetic  (TEM)  sensors  are  most  commonly  used  for 
detection  of  buried  metallic  targets.  These  instruments  actively  transmit  a  time-varying 
primary  magnetic  field  which  illuminates  the  earth.  The  variation  of  the  primary  field 
induces  currents  in  buried  targets  and  these  currents  in  turn  produce  a  secondary  field 
which  can  be  measured  by  a  receiver  at  the  surface  (figure  1) 


Figure  1.  Electromagnetic  induction  (EMI)  survey.  Eddy  currents  are  in¬ 
duced  in  a  buried  target  by  a  time- varying  primary  field.  Decaying  secondary 
fields  radiated  by  the  target  are  then  measured  by  a  receiver  at  the  surface. 

In  a  detection  mode  survey,  a  TEM  sensor  passes  over  an  area  in  nominally  straight, 
parallel  lines,  with  line  spacing  and  instrument  height  dictated  by  instrument  geometry  and 
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detection  considerations.  Subsequent  cued  interrogations  may  revisit  previously-identified 
targets  and  acquire  finely  spaced,  high  signal-to-noise  ratio  (SNR)  data  in  a  small  area 
about  the  target.  Recently  developed  systems  for  cued- interrogation  illuminate  the  tar¬ 
get  with  multiple  transmitters  and  receivers  (a  “multistatic”  configuration)  from  a  single 
observation  location  and  thereby  avoid  the  requirement  for  accurate  positioning  of  mov¬ 
ing  sensors.  Table  1  compares  the  industry-standard  EM-61  sensor  with  newer  multistatic 
systems  developed  specifically  for  UXO  detection  and  classification. 

2.2.  Feature  extraction.  Once  target  anomalies  have  been  identified  in  the  observed  geo¬ 
physical  data,  we  can  characterize  each  anomaly  by  estimating  features  which  will  sub¬ 
sequently  allow  a  classification  algorithm  to  discern  targets  of  interest  (TOI)  from  non- 
hazardous  clutter  (non-TOI).  These  features  may  be  directly  related  to  the  observed  data 
(e.g.  anomaly  amplitude  at  the  first  time  channel),  or  they  may  be  the  parameters  of  a  phys¬ 
ical  model.  Advanced  classification  relies  upon  physics-based  modeling,  with  the  observed 
magnetic  field  B  (t)  radiated  by  a  buried  target  usually  represented  as  a  time- varying  dipole 

(1)  ^(r,t)  =  h!i(3(p(().r)f-p(t)) 

with  r  =  rr  the  separation  between  target  and  observation  location,  and  p(f)  =  p(t)p(t)  a 
time-varying  dipole  moment 

(2)  p (t)  =  -P (t)  •  B0. 

ho 

The  induced  dipole  is  the  projection  of  the  primary  field  B0  onto  the  target’s  polarizability 
tensor  PR)  (Bell  et  ah,  2001).  Here  the  elements  of  the  polarizability  tensor  (PtJ(t))  represent 
the  convolution  of  the  target’s  B-field  impulse  response  (PR))  with  the  transmitter  waveform 
i(t)  (Wait,  1982) 

OO 

(3)  P,j(t)  =  FtJ  -  tW)dH. 

— OO 

The  polarizability  tensor  is  assumed  to  be  symmetric  and  positive  definite  and  so  can  be 
decomposed  as 

(4)  PR)  =  ArLR)A 

with  A  an  orthogonal  matrix  which  rotates  the  coordinate  system  from  geographic  coor¬ 
dinates  to  a  local,  body  centered  coordinate  system.  The  diagonal  eigenvalue  matrix  LR) 
contains  the  principal  polarizabilities  L;R)  (?'  =  1,  2, 3),  which  are  assumed  to  be  independent 
of  target  orientation  and  location. 

From  a  set  of  observations  of  the  electromagnetic  field,  the  inverse  problem  is  then  to 
find  the  set  of  model  parameters  (location,  orientation,  and  polarizabilities)  that  best  fits 
the  data.  The  model  vector  m  can  be  estimated  by  minimizing  a  norm  (e.g.  least  squares) 
quantifying  the  misfit  between  observed  (do6s)  and  predicted  (ydpred)  data.  For  the  TEM 
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Sensor 


Geometry 


Channels 

4  channels 


0.5 


0.040.1  1  10  50 

Time  (ms) 


32  channels 


Table  1.  Electromagnetic  sensors  used  for  UXO  classification.  In  geometry 
plots  transmitters  and  receivers  are  red  and  black,  respectively. 
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dipole  model,  the  least  squares  estimate  must  generally  be  obtained  iteratively,  owing  to  the 
nonlinear  relationship  between  model  parameters  and  predicted  data  in  equation  1.  However, 
if  the  location  of  the  target  (r)  is  assumed  known,  then  the  forward  modeling  becomes  linear, 
so  that 

(5)  dpred  =  Gm 

with  G,  the  forward  modeling  matrix,  implicitly  dependent  on  target  location.  The  least 
squares  model  estimate  is  then  given  by 

m  =  (GTG)_1GTdo6s 

=  Gtdo6S 

with 

(7)  G*  =  (GtG)“1Gt 

denoting  the  pseudo-inverse.  When  inverting  observed  field  data  for  a  sensor  with  tri-axial 
transmit  and  receive  coils  (e.g.  MetalMapper) ,  we  express  G  as 

BsB; 

Eg  B:p  +  B;B; 

ByBy 

BysBp  +  BlBv 

bzsb; 

with  Bp  the  primary  field  at  the  target  and  Bs  the  secondary  field  at  the  receiver,  with 
all  fields  implicitly  dependent  upon  target  (r)  and  sensor  location.  Superscripts  denote  the 
x,  y,  z  components  of  the  respective  fields.  In  this  formulation,  the  model  vector  at  location 
r  is  parameterized  in  terms  of  the  six  unique  elements  of  the  polarizability  tensor  P 

(9)  m  [-fxxj  BXyi  Pxz,  Pyyi  Pyzi  Pzz\ 

In  practice,  the  vector  m  is  estimated  at  each  time  channel  in  a  sequential,  or  two  stage, 
inversion  strategy  (Song  et  al.,  2011).  As  illustrated  in  figure  2,  we  first  solve  a  nonlinear 
inverse  problem  for  target  location.  We  then  solve  a  linear  problem  for  the  polarizability 
tensor  elements  (equation  9)  at  our  fixed  location  estimate.  Decoupling  the  location  and  po¬ 
larizability  parameters  in  this  way  allows  for  efficient  parallel  solution  of  the  linear  problem 
at  all  time  channels.  Target  orientation  and  principal  polarizabilities  are  subsequently  esti¬ 
mated  via  joint  diagonalization  (Cardoso,  1996).  This  algorithm  returns  a  single  eigenvector 
matrix  A  for  all  channels,  corresponding  to  a  fixed  target  orientation.  The  eigenvalues  at 
each  time  channel  are  then  an  estimate  of  principal  polarizabilities. 
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Figure  2.  Sequential  inversion  approach  for  estimation  of  dipole  model  pa¬ 
rameters.  We  first  estimate  target  location  r;  the  predicted  data  in  this  case 
are  a  nonlinear  functional  F( r).  We  then  estimate  target  orientation  and 
polarizabilities.  At  a  fixed  location  and  orientation,  the  predicted  data  are 
related  to  the  model  via  a  linear  forward  operator  G. 

The  model  vector  at  each  time  channel  can  also  be  comprised  of  the  principal  polarizabil¬ 
ities 

(10)  m  =  [Lu  L2,  L3] . 

In  this  case  the  matrix  G  depends  implicitly  upon  both  target  location  and  orientation.  As 
will  be  discussed  in  section  4,  this  formulation  of  the  model  vector  and  forward  operator  is 
useful  when  propagating  noise  on  the  data  to  predicted  distributions  of  polarizabilities. 

2.3.  Classification.  To  rank  detected  targets  for  digging,  we  use  the  information  in  our 
observed  geophysical  data.  Features  of  the  observed  data,  estimated  without  resorting  to 
inversion  with  a  physics-based  model,  can  sometimes  suffice  as  criteria  to  discriminate  be¬ 
tween  ordnance  and  non-ordnance  targets  (e.g.  Williams  et  al.  (2007)).  However,  because 
dipole  model  parameters  can  be  related  to  intrinsic  properties  such  as  target  size  and  shape, 
features  derived  from  the  estimated  parameters  (such  as  the  total  polarizability  in  equation 
3)  are  often  more  reliable  for  discriminating  between  TOI  and  non-TOI. 

Classification  with  TEM  data  is  often  performed  by  comparing  estimated  polarizabilities 
with  library  responses  and  then  ranking  a  target  based  on  some  measure  of  closeness  between 
observed  and  expected  responses.  Care  must  be  taken  here  to  use  parameters  which  can 
be  reliably  estimated:  late  time  polarizabilities  are  more  susceptible  to  noise  and  poor 
estimates  may  unduly  affect  the  discrimination  decision.  Pasion  et  al.  (2007)  solve  this 
problem  with  a  fingerprinting  algorithm  that  inverts  for  target  location  and  orientation  while 
holding  principal  polarizabilities  fixed  at  their  library  values.  Reducing  the  model’s  degrees 
of  freedom  in  this  way  makes  the  inversion  less  susceptible  to  fitting  the  noise.  Targets  are 
then  dug  based  upon  the  proposed  library  item  which  produces  the  best  fit  to  the  observed 
data. 

Statistical  classification  algorithms  applied  to  UXO  discrimination  try  to  learn  a  decision 
rule  from  a  sample  of  labeled  targets  for  which  ground  truth  is  known  (the  training  data). 
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One  approach  to  formulating  the  decision  rule  is  to  fit  some  assumed  parametric  distribu¬ 
tions  to  each  class  of  targets  in  the  training  data,  and  to  then  assign  an  unlabeled  (test) 
target  to  the  class  distribution  which  is  most  likely.  The  class  distributions  are  defined  in  a 
multidimensional  feature  space  spanned  by  some  subset  of  estimated  model  parameters,  or 
transformations  thereof.  The  total  polarizability 

3 

(ii)  =  *;) 

2=1 

can  be  a  useful  classification  feature  when  secondary  and  tertiary  (L2,  L:>)  polarizabilities 
are  poorly  constrained.  Alternatively,  a  two-dimensional  space  spanned  by  the  amplitude 
and  decay  of  the  polarizabilities 


has  been  successfully  used  to  train  statistical  classifiers  for  a  number  of  ESTCP  field  demon¬ 
strations  (e.g.  Billings  et  al.  (2010)).  These  parameters  are  useful  because,  to  first  order, 
a  conductor  can  be  modeled  as  a  simple  LR  (circuit  with  inductor  and  resistor  in  series) 
loop  which  is  inductively  coupled  to  transmitters  and  receivers  on  the  surface.  The  current 
response  of  this  loop  is  a  decaying  exponential  which  is  fully  described  by  an  amplitude  and 
time  constant  (West  and  Macnae,  1991).  Polarizabilities  estimated  from  multistatic  data  are 
typically  well-constrained  relative  to  monostatic  detection  data,  and  so  statistical  classifiers 
for  cued  data  sets  are  usually  trained  on  the  (log-transformed)  principal  polarizabilities,  or 
total  polarizabilities. 

The  output  of  any  automated  classification  algorithm  is  a  decision  statistic,  or  score, 
that  is  used  to  rank  detected  targets  from  likely  TOI  to  likely  non-TOI.  For  example,  a 
library  classifier  uses  the  misfit  of  estimated  polarizabilities  with  library  polarizabilities  as 
the  decision  statistic.  As  shown  in  Figure  3,  the  receiver  operating  characteristic  (ROC)  is 
then  a  plot  of  the  true  positive  fraction  ( TPF )  versus  the  false  positive  fraction  (. FPF ), 
which  are  defined  as  the  cumulative  score  distributions  for  TOI  and  non-TOI 


(13) 


TPF(x ) 


J  p(y\TOI)dy 


FPF(x) 


J  p(y\non  —  TOI)dy. 


In  the  context  of  munitions  management,  the  false  alarm  rate  (FAR)  at  which  all  ordnance 
are  detected  on  the  ROC  (i.e.  TPF  =  1)  is  the  crucial  metric  by  which  site  managers 
evaluate  the  efficacy  of  remediation  efforts.  An  advanced  technique  that  results  in  good 
initial  detection  of  TOI  but  fails  to  find  outlying  TOI  until  late  in  the  dig  order  will  be 
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Figure  3.  Left:  the  receiver  operating  characteristic  (ROC)  shows  the  true 
positive  fraction  (TPF)  as  a  function  of  the  false  positive  fraction  (FPF). 

Right:  the  ROC  at  point  x  can  be  modeled  as  the  integral  of  the  distributions 
of  TOI  and  non-TOI  (true  and  false  positives,  respectively)  with  respect  to 
the  decision  statistic. 

judged  unsuccessful.  This  can  occur  with  statistical  classifiers  that  overfit  the  training  data 
and  so  fail  to  generalize  to  noisy  test  TOI. 

3.  Predicting  target  response  thresholds 

Targets  identified  in  a  detection  mode  survey  must  be  selected  for  subsequent  cued  in¬ 
terrogation  using  objective  and  reproducible  criteria.  The  Geophysical  System  Verification 
(GSV)  approach  developed  by  ESTCP  (Nelson  et  ah,  2010)  advocates  setting  a  target  pick¬ 
ing  threshold  based  upon  the  lowest  predicted  amplitude  expected  for  targets  of  interest 
at  a  specified  clearance  depth.  In  the  case  of  the  EM-61  cart,  the  worst  case  scenario  is  a 
horizontal  item,  oriented  perpendicular  to  the  sensor  track.  This  geometry  decouples  the 
primary  field  from  the  axial  (Li)  response  of  the  target,  so  that  the  dipole  response  is  dom¬ 
inated  by  the  sum  of  the  smaller,  transverse  polarizabilities.  The  GSV  response  calculator 
provides  functionality  for  predicting  the  best  and  worst  case  target  response  amplitudes  for 
an  array  of  EM-61  sensors  in  an  arbitrary  geometry. 

When  we  apply  this  approach  to  munitions  response  projects,  we  sometimes  find  that 
the  threshold  defined  by  the  worst  case  scenario  results  in  an  undesirably  large  number  of 
target  picks.  This  has  motivated  us  to  extend  the  GSV  response  calculator  functionality  to 
consider  the  statistics  of  the  target  response  across  the  swath  of  the  array. 

Figure  4  illustrates  this  extended  analysis  for  a  three  coil  EM-61  array  interrogating  a 
37  mm  projectile  at  30  cm  depth.  Within  the  swath  of  the  array,  the  minimum  predicted 
anomaly  amplitude  occurs  for  orientations  of  the  target  across  the  sensor  track.  The  global 
minimum  of  the  response  (~  8.5  mV  at  channel  3)  is  near  the  cross-track  edges  of  the 
center  transmitter  coil.  This  minimum  is  the  consequence  of  the  fall-off  in  received  EMF 
as  the  secondary  dipole  moves  away  from  the  center  of  each  receiver  coil.  In  this  particular 
example,  the  minimum  response  is  in  a  cross  track  orientation. 

For  asymmetric  array  geometries  we  have  found  that  the  minimum  response  is  still  at 
a  horizontal  dip,  but  the  target  azimuth  can  be  at  an  intermediate  angle  to  cross  and 
along  track  directions.  This  is  illustrated  in  figure  5  for  an  asymmetric  three  coil  array. 
When  constructing  the  worst  case  curve  we  therefore  model  over  a  full  360  degree  range 
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of  target  azimuths  to  obtain  the  minimum  predicted  response.  Comparing  the  symmetric 
and  asymmetric  array  geometries  in  figures  4  and  5,  we  remark  that  the  minimum  response 
across  the  swath  is  significantly  elevated  for  the  asymmetric  array.  Narrowing  the  swath  will 
decrease  the  production  rate,  but  will  improve  target  detection. 


(a)  Along  track  orientation 


(b)  Across  track  orientation 


Cross  track  (m) 


(c)  Minimum  response 


Figure  4.  Anomaly  amplitude  analysis  for  a  horizontal  37  mm  projectile  at 
30  cm  depth  illuminated  by  a  three  coil  EM-61  array,  (a)  Maximum  anomaly 
amplitude  (mV)  at  channel  3  as  a  function  of  target  location,  for  target  ori¬ 
ented  along  track.  White  lines  indicate  transmitters,  (b)  As  in  (a),  but  with 
horizontal  target  oriented  across  track,  (c)  Minimum  of  responses  in  (a)  and 
(b)  across  the  swath  of  the  sensor  array. 

Returning  now  to  the  symmetric  three  coil  array,  we  note  that  the  8.5  mV  minimum  in 
figure  4(c)  will  undoubtedly  produce  a  very  large  number  of  picks  that  must  be  excavated  or 
interrogated  with  a  cued  sensor.  For  example,  at  the  Jeep  and  Demolition  Range  (JDR)  at 
the  Former  Lowry  Bombing  and  Gunnery  Range  (FLBGR),  we  found  that  a  30  cm  clearance 
depth  for  37  mm  projectiles,  corresponding  to  an  8.5  mV  threshold,  required  us  to  dig  a 
prohibitively  large  number  of  picks. 

However,  we  have  found  that  by  considering  the  statistics  of  the  target  response  we  can 
significantly  reduce  the  number  of  target  picks.  Assuming  a  uniform  spatial  distribution  of 
horizontal  37  mm  projectiles  throughout  the  sensor  footprint  with  a  uniform  distribution  of 
target  azimuths,  we  obtain  the  cumulative  distribution  of  response  amplitudes  in  figure  6. 
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(a)  Along  track  orientation 


-1  0  1 

Cross  track  (m) 


(b)  Across  track  orientation 


-1  0  1 

Cross  track  (m) 


(c)  Minimum  response 


Figure  5.  Anomaly  amplitude  analysis  for  an  asymmetric  three  coil  EM- 
61  array  Note  in  this  case  that  the  worst  case  amplitude  does  not  always 
correspond  to  the  cross  track  orientation. 


At  0.99  and  0.95  confidence  levels,  we  can  raise  the  detection  threshold  in  figure  4  to  ap¬ 
proximately  11  and  15  mV,  respectively.  In  (c)  we  plot  the  proportion  of  37  mm  targets  that 
fall  below  these  thresholds  as  a  function  of  cross  track  location.  For  the  11  mV  threshold 
(0.99  confidence  level),  the  non-zero  regions  of  targets  below  the  threshold  are  restricted  to 
lobes  near  the  outer  edges  of  the  center  transmitter.  Increasing  the  threshold  to  15  mV  (0.95 
confidence  level)  significantly  increases  the  chances  that  we  will  miss  a  37  mm  target  within 
the  sensor  swath. 

As  shown  in  figure  7,  increasing  the  threshold  to  11  mV  at  JDR  reduced  the  number 
of  target  picks  by  30%,  while  still  ensuring  a  very  high  confidence  of  finding  all  37  mm 
projectiles.  At  the  95%  confidence  threshold  we  have  a  60%  reduction  in  target  picks. 
However,  the  5%  probability  of  missing  a  horizontal  37  mm  may  not  be  an  acceptable  risk. 
At  sites  where  the  worst  case  threshold  produces  a  prohibitive  number  of  picks,  this  type 
of  statistical  analysis  can  provide  objective  justification  for  slightly  raising  the  detection 
threshold. 
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(a)  (b) 


Minimum  response  (mV)  Minimum  response  (mV) 


Figure  6.  (a)  Cumulative  distribution  of  response  amplitudes  for  horizontal 
37  mm  projectiles  at  30  cm  depth.  Targets  are  distributed  uniformly  with 
respect  to  cross  track  location  and  azimuth,  (b)  Lower  tail  of  the  cumulative 
distribution.  Markers  indicate  0.99  and  0.95  confidence  bounds,  (c)  Propor¬ 
tion  of  horizontal  37  mm  targets  with  response  below  threshold  specified  by 
confidence  level  1  —  a,  as  a  function  of  cross  track  location. 

The  distribution  in  figure  6  is  generated  for  horizontal  37  mm  targets  at  arbitrary  az¬ 
imuths.  Allowing  for  an  arbitrary  target  azimuth  and  dip  generates  a  somewhat  less  con¬ 
servative  distribution  of  response  amplitudes.  Figure  8  shows  the  empirical  probability  and 
cumulative  density  functions  for  simulations  assuming  a  uniform  angular  distribution  of  tar¬ 
gets.  These  results  are  again  generated  for  37  mm  projectiles  at  30  cm  depth,  with  the 
symmetric  EM-61  array  of  figure  4. 

In  future  work  we  will  extended  this  type  of  GSV  functionality  to  newer  dynamic  plat¬ 
forms,  including  the  MetalMapper  and  TEMTADS2x2.  If  requested  by  the  program  office, 
we  will  work  to  develop  a  software  tool,  similar  to  the  existing  GSV  calculator,  that  makes 
this  analysis  widely  available  to  the  munitions  response  industry. 
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Figure  7.  Cumulative  distribution  of  response  amplitudes  at  JDR.  Markers 
indicate  1  —  a  confidence  thresholds  for  detection  of  37  mm  targets. 
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Figure  8.  Distribution  of  response  amplitudes  for  37  mm  targets  at  30  cm 
depth,  distributed  uniformly  in  azimuth  and  dip.  (a)  Empirical  distribution. 
Vertical  lines  are  1  —  a  confidence  level  cut-off.  (b)  Cumulative  distribution, 
(c)  Proportion  of  37  mm  targets  with  response  below  threshold  specified  by 
1  —  a  confidence  level,  as  a  function  of  cross  track  location. 
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4.  Performance  prediction 

Excellent  results  have  been  obtained  throughout  the  ESTCP  demonstration  program  with 
the  sensor  and  data  processing  technologies  described  above:  advanced  classification  con¬ 
sistently  outperforms  conventional  processing  of  industry-standard  detection  data  (Billings 
et  al.  (2010),  Steinhurst  et  al.  (2010),  Prouty  et  al.  (2011),  Shubitidze  et  al.  (2011)).  The 
significant  reductions  in  the  false  alarm  rate  obtained  with  advanced  classification  translate 
to  substantial  cost-savings  during  site  remediation.  Motivated  by  these  successes,  classifica¬ 
tion  methods  have  increasingly  been  transitioned  from  researchers  to  industry.  For  example, 
at  the  2011  Pole  Mountain  demonstration,  production  geophysicists  working  with  UBC- 
GIF  software  to  process  MetalMapper  data  obtained  classification  performance  comparable 
to  that  of  “expert”  analysts  (Pasion,  2012).  These  initial  results  indicate  that  advanced 
classification  can  be  fully  transitioned  to  the  munitions  response  industry. 

To  further  this  transition  process,  we  have  developed  methods  for  efficient  prediction  of 
classification  performance.  The  goal  of  this  work  is  to  provide  tools  for  site  managers  and 
their  teams  to  reliably  and  objectively  predict  the  benefits  of  advanced  classification  for  a 
particular  remediation  problem.  In  particular,  we  address  two  questions: 

(1)  Will  multi-static  sensor  data  reliably  discriminate  between  targets  of  interest  and 
representative  clutter  items  encountered  at  the  site? 

(2)  Is  there  a  significant  improvement  in  expected  classification  performance  -  and  hence 
a  reduction  in  expected  costs  -  relative  to  classification  with  a  monostatic  sensor  such 
as  the  EM-61?  For  a  relatively  simple  classification  task  (e.g.  identification  of  4.2” 
mortars  at  Camp  Sibert,  Billings  et  al.  (2010)),  the  additional  survey  and  process¬ 
ing  costs  associated  with  advanced  classification  may  ultimately  provide  negligible 
improvements  in  classification. 

Throughout  this  report  we  will  exploit  available  data  sets  from  ESTCP  demonstrations  to 
develop  and  validate  our  approach  to  performance  prediction.  A  summary  of  these  data  sets 
is  provided  in  Appendix  A. 

As  discussed  in  the  previous  section,  survey  design  for  munitions  response  projects  is 
largely  restricted  to  detection  considerations.  Here  we  shift  the  emphasis  in  survey  design 
from  detection  to  classification.  We  consider  how  site-specific  inputs  (sensor  platform,  tar¬ 
get  type,  etc.)  ultimately  affect  classification  performance,  as  quantified  by  the  ROC.  Two 
routes  can  be  taken  to  map  from  these  inputs  to  a  predicted  ROC  curve.  First,  numerical 
simulation  of  synthetic  observed  sensor  data  can  be  used  to  quantify  the  variability  of  es¬ 
timated  polarizabilities,  and  hence  to  predict  the  score  distributions  of  TOI  and  non-TOI 
for  a  given  noise  model.  This  approach  is  numerically  intensive  and  is  not  suited  for  an 
efficient  decision  support  tool.  We  therefore  pursue  the  second  option:  we  develop  statis¬ 
tical  mappings  between  environmental  inputs  and  estimated  polarizabilities,  as  depicted  in 
figure  9. 
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Figure  9.  Methods  for  polarizability  prediction.  Most  studies  use  inversion 
of  synthetic  data  with  independent,  Gaussian  noise  to  predict  polarizability 
distributions.  In  the  next  section,  we  develop  statistical  mappings  that  charac¬ 
terize  correlated,  non-Gaussian  noise  with  a  noise  covariance,  then  propagate 
noise  to  the  model  and  integrate  over  specified  spatial  and  angular  distribu¬ 
tions  of  targets. 

4.1.  Performance  prediction  for  cued  interrogation. 

4.1.1.  Noise  estimation.  We  begin  by  quantifying  the  noise  on  the  observed  data  for  a  sta¬ 
tionary  cued  sensor  such  as  the  MetalMapper.  The  problem  is  somewhat  simplified  by  the 
fact  that  soundings  are  acquired  at  a  single  location;  we  need  not  account  for  the  errors  in 
measurement  of  sensor  location  that  complicate  analysis  of  dynamically-collected  data.  The 
errors  on  measured  cued  data  are  then  a  combination  of  sources,  including  baseline  sensor 
noise  and  EMI  responses  that  are  unaccounted  for  in  an  inversion  (e.g.  transient  cultural 
noise,  magnetic  soils  or  neighboring  targets). 

Background  noise  is  characterized  via  repeated  measurements  taken  at  an  electromagneti- 
cally  quiet  location  at  the  site.  These  measurements  quantify  sensor  noise  and  soil  response. 
Significant  variability  in  sensor  noise  was  observed  with  two  MetalMapper  systems  used 
for  the  2010  Camp  Butner  demonstration.  The  older,  prototype  system  had  substantially 
higher  noise  levels  than  the  newer  system  manufactured  by  Geometries.  This  difference  was 
likely  due  to  improvements  in  electronics  in  the  newer  system.  Subsequent  deployments  of 
Geometries  MetalMappers  at  Pole  Mountain  and  Camp  Beale  have  produced  consistently 
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low  noise  levels  associated  with  the  sensor  platform  and  electronics.  Figure  10  compares 
standard  deviations  estimated  from  background  measurements  taken  at  Camp  Butner  and 
Pole  Mountain. 


Figure  10.  Background  noise  standard  deviations  for  MetalMapper  sensor 
at  Camp  Butner  and  Pole  Mountain. 


When  inverting  field  data,  we  assume  that  the  observed  data  are  a  superposition  of  the 
true  data  predicted  by  one  or  more  dipole  sources,  and  a  noise  term  e 

(14)  dobs  =  dtrue  +  e 


with  e  ~  A/"(0,  a,)  a  zero  mean,  normally-distributed  random  variable.  The  estimated  stan¬ 
dard  deviations  weight  the  contributions  of  each  channel  to  the  data  misfit 

'd°b‘  -  2 


(15) 


fa  =  52 
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In  the  sequential  inversion  approach  depicted  in  figure  2,  this  weighting  is  less  critical  since 
each  channel  is  inverted  separately  at  a  fixed  location.  Channel  weighting  is  more  important 
when  all  channels  are  inverted  at  once,  as  in  the  fingerprinting  approach  described  in  Pasion 
(2007). 

Of  course,  the  estimated  standard  deviations  derived  from  background  measurements  do 
not  account  for  the  various  non-stationary  sources  of  noise  that  are  inevitably  present  in  field 
data.  Figure  11  shows  one  such  source:  an  elevated  magnetic  soil  response  encountered  in 
part  of  the  survey  area  at  Camp  Butner.  Non-dipolar  fields  that  cannot  be  fit  in  an  inversion 
can  also  be  regarded  as  additive  noise  on  the  observed  data.  Given  the  complex  dependence 
of  noise  upon  sensor  and  environment,  we  do  not  attempt  to  predict  individual  noise  sources 
and  their  subsequent  effect  on  estimated  polarizabilities.  Rather,  we  use  the  distributions  of 
polarizabilities  recovered  from  ESTCP  field  data  sets  to  estimate  the  effective  noise  across 
each  site. 

Figure  12  compares  polarizabilities  recovered  for  small  industry  standard  objects  (ISOs) 
at  Pole  Mountain  and  Camp  Beale.  It  is  evident  from  the  distributions  of  polarizabilities 
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Figure  11.  Components  of  the  observed  secondary  field  for  MetalMapper 
dynamic  data  acquired  at  Camp  Butner.  Regions  of  elevated  magnetic  soil 
response  are  delineated  in  the  z-component  receiver  data.  Horizontal  compo¬ 
nents  (x  and  y)  of  the  data  are  less  sensitive  to  magnetic  soils  when  excitation 
is  from  a  z-component  transmitter  only. 


that  the  effective  noise  at  Camp  Beale  is  larger  than  at  Pole  Mountain.  Assuming  perfect 
recovery  of  target  location,  each  set  of  estimated  polarizabilities  in  Figure  12  is  a  linear 
transformation  of  the  observed  data  acquired  over  an  ISO  target  (see  equation  6).  The 
covariance  of  the  jth  polarizability  model  at  the  ith  time  channel  can  then  be  expressed  as 

(16)  covfmy)  =  (G*)rcov(df*)G] 
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Figure  12.  Estimated  ISO  polarizabilities  at  Camp  Beale  and  Pole  Mountain. 
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with  the  pseudo-inverse  Gj  given  by  equation  7.  For  uncorrelated  errors  with  constant 
standard  deviation  cq  at  each  time  channel,  the  data  covariance  is  cov(df)S)  =  of  I  and  the 
model  covariance  simplifies  to 

(17)  cov(niy)  =  oftGjG,-)"1 

As  discussed  in  section  2.2,  the  forward  modeling  operator  G  can  be  computed  for  a  fixed 
target  location  and  orientation,  in  which  case  the  model  vector  at  each  time  channel  is 
composed  of  the  principal  polarizabilities  at  that  channel  (equation  10). 

The  covariance  of  the  principal  polarizabilities  in  equation  17  represents  the  covariance  of 
the  model  parameters  that  would  be  estimated  from  repeated  noise  realizations  with  stan¬ 
dard  deviation  cq  for  a  single  target  at  a  fixed  location  and  orientation.  The  polarizabilities 
estimated  from  field  data  are  instead  a  set  of  model  estimates  with  independent  realizations 
of  noise  but  at  varying  locations  and  orientations.  The  following  procedure  computes  an  es¬ 
timator  of  the  effective  noise  standard  deviation  at  a  site  from  the  set  of  recovered  principal 
polarizabilities: 


Algorithm  1  Estimator  of  of  the  data  variance. 

(1)  Estimate  the  (3  x  3)  sample  covariance  matrix  using  the  recovered  principal  polariz¬ 
abilities  at  ith  time  channel,  for  all  targets  within  a  given  class. 

(2)  Form  a  vector  h,  comprised  of  the  diagonal  elements  of  the  estimated  model  covari¬ 
ance  (i.e.  the  variances  of  the  polarizabilities) 

(3)  At  each  target’s  estimated  location  and  orientation,  compute  IT,  =  (GjGj)-1.  This 
is  the  model  covariance  for  the  jth  target  (equation  17)  without  scaling  by  the  (un¬ 
known)  data  variance  of 

(4)  Compute  the  mean  model  covariance  matrix 

1  M 

<18>  r=M £r* 

3= 1 

over  all  M  targets  in  the  class.  T  is  an  average  of  symmetric  and  positive-definite 
covariance  matrices  and  so  is  itself  a  valid  covariance  matrix. 

(5)  Form  the  vector  g  from  the  diagonal  elements  of  f . 

(6)  Finally,  compute  a  least  squares  estimate  of  the  data  variance  at  the  ith  time  channel 

(19)  a?  =  (gTg)-1gTh, 

To  ensure  that  the  inverse  problem  in  equation  19  is  well-posed,  we  assume  a  diago¬ 
nal  data  covariance  and  neglect  estimation  of  correlations  quantified  by  off-diagonal 
terms. 


Figure  13  shows  a  synthetic  example  of  noise  estimation  for  MetalMapper  data.  The  syn¬ 
thetic  observed  data  are  contaminated  with  additive  zero  mean  Gaussian  noise  with  standard 
deviations  derived  from  Camp  Butner  background  measurements.  We  simulate  ISOs  dis¬ 
tributed  uniformly  in  both  location  and  orientation,  with  a  maximum  target  depth  of  0.6 
m.  The  recovered  distribution  of  polarizabilities,  estimated  from  nonlinear  inversion  of  the 
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synthetic  data,  is  then  used  to  compute  the  data  variance  estimator  using  the  above  proce¬ 
dure.  There  is  good  agreement  between  the  actual  and  estimated  data  standard  deviations, 
suggesting  that  this  approach  can  provide  a  reasonable  estimate  of  the  effective  noise  on  the 
data. 
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Figure  13.  Synthetic  example  of  estimation  of  effective  noise  standard  devi¬ 
ation.  Left:  ISO  polarizabilities  estimated  from  synthetic  MetalMapper  data 
contaminated  with  Gaussian  noise.  The  kink  in  the  true  secondary  polariz¬ 
abilities  at  late  times  is  a  consequence  of  noise  on  the  particular  measurements 
used  to  extract  reference  polarizabilities.  Smoothing  could  be  applied  to  re¬ 
move  this  effect.  Right:  True  standard  deviation  of  noise  and  estimator  d\ 
recovered  from  polarizability  distributions. 


We  now  apply  this  analysis  to  real  MetalMapper  data.  The  analysis  is  restricted  to  ESTOP 
demonstration  sites  with  seeded  ISOs  (at  the  time  of  this  report):  Pole  Mountain  and  Camp 
Beale.  ISO  items  are  preferred  for  this  analysis  because  of  their  physical  uniformity  and 
consequent  consistency  of  their  polarizabilities  under  low  noise  conditions.  ISOs  can  also  be 
seeded  uniformly  across  the  site  to  ensure  that  variability  in  site  conditions  is  adequately 
sampled. 

Figure  14  shows  noise  estimates  for  Camp  Beale  and  Pole  Mountain  derived  from  the 
distributions  of  estimated  ISO  polarizabilities  in  figure  12.  As  expected  from  the  observed 
distributions  of  polarizabilities,  the  estimated  standard  deviation  of  the  noise  at  Camp  Beale 
is  larger  than  at  Pole  Mountain,  but  only  slightly.  Munkholm  and  Auken  (1996)  showed  that 
log-gated  and  stacked  white  noise  produces  errors  on  the  TDEM  response  with  a  standard 
deviation  exhibiting  a  t~1'2  decay.  For  these  data  sets,  the  estimated  noise  standard  devia¬ 
tions  decay  considerably  faster  than  t-1. 

As  a  check  on  this  procedure,  we  can  now  contaminate  synthetic  data  with  our  effective 
noise  and  verify  that  the  distributions  of  recovered  polarizabilities  are  consistent  with  the 
actual  distributions  in  figure  12.  Figure  15  compares  the  synthetic  estimates  recovered  with 
the  inferred  Pole  Mountain  noise  model  and  the  actual  estimated  polarizabilities  for  ISOs 
at  this  site.  The  synthetic  polarizabilities  have  a  much  higher  variance,  in  particular  for  the 
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secondary  and  tertiary  polarizability  estimates,  indicating  that  we  have  overestimated  the 
noise.  A  similar  result  is  obtained  for  simulations  with  our  estimated  Camp  Beale  noise. 

Why  does  our  noise  estimation  procedure  fail  in  this  case?  We  have  made  several  as¬ 
sumptions  in  this  analysis:  the  noise  is  independent,  normally  distributed,  and  constant  for 
all  data  at  each  time  channel.  As  we  will  now  show,  a  closer  examination  of  the  polariz¬ 
abilities  estimated  from  field  data  indicates  that  the  noise  on  the  data  is  correlated,  only 
approzimately  normally-distributed,  and  is  not  constant  for  all  data  measured  at  a  single 
time  channel. 

The  correlated  nature  of  the  noise  is  easy  to  see  in  figure  12:  correlated  shifts  of  the  esti¬ 
mated  polarizabilities,  especially  obvious  at  Camp  Beale,  could  not  result  from  independent 
random  noise.  In  the  case  of  Camp  Beale,  a  variable  magnetic  soil  response  affects  our 
estimation  of  target  depth  and  so  results  in  constant  shifts  (in  log  space)  of  the  estimated 
polarizabilities. 

We  turn  now  to  the  distribution  of  the  noise.  Typically,  in  synthetic  simulation  work 
we  assume  Gaussian  noise  at  each  time  channel.  This  is  motivated  by  simplicity,  and  the 
observation  that  multiple  sources  of  additive  noise  will  tend  to  a  normal  distribution. 

In  previous  work,  we  also  carried  out  simulations  that  treated  the  data  residuals  -  the 
difference  between  observed  and  predicted  data  -  as  realizations  of  noise  (Carin  et  ah,  2012). 
The  various  realizations  of  noise  across  a  site  were  then  added  to  forward  modeled  data  for 
a  given  target  to  synthetically  seed  that  target  at  the  site.  This  allowed  us  to  study  the 
variability  of  estimated  target  features  expected  at  the  site. 

A  drawback  of  the  synthetic  seeding  method  is  that  the  residuals  in  an  inversion  in  fact 
represent  the  component  of  the  noise  that  is  not  fit  in  an  inversion.  We  are  not  accounting 
for  the  component  of  the  noise  that  causes  random  deviations  of  the  estimated  model  and 
contributes  to  the  overall  variability  of  the  features  within  a  target  class.  Here  we  try  to 


Figure  14.  Estimates  of  effective  noise  standard  deviations  for  MetalMapper 
Beale  Parsons  and  Pole  Mountain  data  sets.  Heavy  dashed  line  is  a  f-1?2 
decay. 
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Figure  15.  Comparison  of  actual  (left)  and  synthetic  estimated  (right)  ISO 
polarizabilities  at  Pole  Mountain.  The  latter  are  generated  by  adding  the 
predicted  data  and  independent  Gaussian  noise  with  standard  deviations  es¬ 
timated  with  algorithm  1. 


isolate  this  effect  of  the  noise  by  fitting  a  smooth,  physically-motivated  function  to  the 
estimated  polarizabilities.  This  provides  a  route  for  estimating  the  component  of  the  noise 
that  contaminates  our  model  estimates. 

The  polarizability  decay  for  an  arbitrary  target  can  be  generally  parameterized  as  a  su¬ 
perposition  of  decaying  exponentials 

N 

(20)  Lsmooth{tj)  =  ^2  ai  exP(— bf/hi)- 

2=1 

With  N  — >■  oo,  this  is  the  analytic  form  of  the  polarizability  response  of  a  spheroidal  tar¬ 
get  (Kaufman,  1994).  Random  fluctuations  of  recovered  polarizabilities  away  from  this 
smoothly  decaying  function  are  then  diagnostic  of  random  noise  on  the  data.  We  fit  the 
above  function  to  the  recovered  polarizabilities,  with  the  summation  truncated  to  N  =  30 
terms  and  decay  constants  logarithmically-spaced  between  10-6  and  10-2  s.  We  solve 
a  linear  optimization  problem  for  the  coefficients  cq  by  minimizing  the  least  squares  differ¬ 
ence  between  estimated  and  smoothed  polarizabilites,  with  the  constraint  a,  >  0  ensuring 
a  monotonically  decreasing  function.  Figure  16  shows  a  typical  fit  obtained  with  this  ap¬ 
proach.  There  is  excellent  agreement  between  smoothed  and  unsmoothed  polarizabilities, 
the  only  obvious  deviation  is  at  late  times  when  the  estimated  secondary  polarizabilities 
approach  the  noise  floor.  Given  estimated  (unsmoothed)  and  smoothed  polarizabilities,  we 
can  generate  an  estimate  of  the  noise  by  predicting  the  data  for  both 
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Figure  16.  Example  fit  (solid  lines)  of  a  sum  of  exponential  decays  to  esti¬ 
mated  polarizabilities  (markers)  for  an  ISO  target  at  Pole  Mountain 


where  the  model  vectors  m  are  the  polarizabilities  at  each  channel.  Our  prediction  for  the 
realization  of  noise  on  the  data  is  then 

(22)  e  des£  d  smooth  GAm. 

From  this  we  see  that  the  predicted  noise  is  that  of  a  dipole  source,  co-located  with  the 
target  and  with  principal  polarizabilities  given  by  the  change  in  the  polarizability  model 
Am.  Figure  17  shows  the  prediction  for  e  for  one  channel  of  MetalMapper  data,  using  the 
ISO  polarizability  estimate  in  figure  16.  While  this  approach  dictates  dipolar  noise,  the 
actual  realization  of  noise  that  produced  the  estimated  (unsmoothed)  polarizabilities  is  not 
necessarily  dipolar  -  there  may  be  another  noise  model  that  can  produce  the  same  result. 
Nonetheless,  we  can  use  the  statistics  of  the  predicted  dipolar  noise  to  understand  the  char¬ 
acter  of  the  noise  that  produced  the  observed  polarizabilities.  Repeating  the  polarizability 
fitting  process  for  all  ISO  items  at  Pole  Mountain,  we  obtain  the  noise  distribution  for  a 
single  channel  in  figure  18.  The  resulting  distribution  can  be  decribed  as  mostly  normal: 
larger  amplitude  outliers  produce  heavier  tails  than  would  be  obtained  with  a  purely  Gauss¬ 
ian  noise  distribution.  This  type  of  heavy-tailed  noise  often  motivates  the  use  of  robust 
norms  that  downweight  the  influence  of  outliers  on  model  estimates.  We  investigated  robust 
fitting  in  an  earlier  project  and  showed  that  these  methods  are  especially  beneficial  when 
processing  dynamic  data  collected  with  monostatic  EM  sensors  (Beran  et  ah,  2011).  In 
contrast,  the  large  number  of  (nearly)  redundant  measurements  acquired  with  multistatic 
sensors  make  these  data  inherently  robust  to  outliers.  Robust  processing  therefore  offers 
negligible  improvements  in  parameter  estimation  for  these  platforms. 

Symmetric,  heavy-tailed  noise  is  often  represented  as  a  mixture  of  two  Gaussian  compo¬ 
nents  with  equal  means  (Marrona  et  ah,  2006).  The  lower  weight,  larger  variance  component 
describes  the  outlying  noise.  In  figure  18  we  find  that  a  single  Gaussian  does  not  adequately 
describe  the  noise,  and  so  we  use  two  techniques  to  fit  univariate,  Gaussian  mixture  models 
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Figure  17.  MetalMapper  data  at  channel  42  (7.9  ms)  predicted  using  un¬ 
smoothed  (top  row)  and  smoothed  (middle  row)  polarizability  estimates  from 
figure  16.  The  difference  between  the  predicted  data  sets  (bottom  row)  is  it¬ 
self  a  dipolar  anomaly  and  is  an  estimate  of  the  noise  realization  on  the  data. 
Markers  indicate  receiver  locations. 


to  the  empirical  distribution  of  the  noise.  Expectation-maximization  (E-M)  is  a  popular, 
iterative  algorithm  for  maximum  likelihood  estimation  of  mixture  model  parameters  (Hastie 
et  ah,  2001).  We  have  also  developed  an  alternative  estimation  method  that  minimizes  the 
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Figure  18.  Fits  to  noise  distribution  derived  from  ISO  target  polarizabili¬ 
ties,  for  channel  42  of  MetalMapper  data  at  Pole  Mountain.  Left:  Empirical 
density  (histogram)  and  predicted  densities  using  a  single  normal  distribution 
(denoted  Normal)  and  a  mixture  of  normal  distributions.  The  mixture  models 
are  estimated  using  expectation-maximization  (E-M)  and  a  novel  method  that 
directly  fits  the  CDF.  Right:  empirical  (black  line)  and  predicted  cumulative 
densities.  Inset  shows  fits  in  the  shoulder  of  the  distribution. 
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least  squares  difference  between  empirical  (observed)  and  predicted  cumulative  distribution 
functions  (CDFs).  By  specifying  particular  points  on  the  CDF  that  we  want  to  fit,  we  place 
more  emphasis  on  the  “shoulders”  of  the  distribution  and  obtain  better  agreement  between 
observed  and  predicted  cumulative  distributions.  Fitting  the  CDF  directly  is  motivated  by 
the  fact  that  the  Kolmogorov-Smirnov  (K-S)  hypothesis  test  of  two  non-parametric  distri¬ 
butions  uses  the  maximum  difference  between  CDFs  as  its  test  statistic  (Press  et  al.,  1992). 
By  this  criterion,  our  method  produces  a  higher  likelihood  estimate  of  mixture  parameters 
than  E-M.  The  CDF  fitting  method  produces  a  good  fit  across  all  time  channels,  indicating 
that  a  two-component  mixture  is  an  accurate  model  for  the  noise. 

Figure  19  shows  the  standard  deviations  and  weights  of  the  normal  mixture  components 
across  all  MetalMapper  channels  at  Pole  Mountain.  The  lowest  noise  variances  are  at  inter¬ 
mediate  channels.  At  late  times  the  observed  decays  hit  the  background  noise  floor,  so  that 
the  noise  standard  deviation  is  a  larger  proportion  of  the  true  datum. 


Figure  19.  Estimated  normal  mixture  model  parameters  at  each  MetalMap¬ 
per  channel.  Left:  Standard  deviations  of  normal  mixture  components,  ex¬ 
pressed  as  a  percentage  of  the  observed  datum,  estimated  via  CDF  fitting 
technique.  Right:  Weighting  of  first,  smaller  variance  component  in  mixture 
model. 


In  this  section  we  have  investigated,  in  some  detail,  methods  for  estimating  the  distribution 
of  additive  noise  that  contaminates  observed  data.  We  conclude  that  a  simple  model  of 
uncorrelated  Gaussian  noise  is  only  roughly  applicable  to  MetalMapper  data.  A  mixture  of 
correlated  Gaussians  is  a  much  better  approximation.  With  these  lessons  in  mind,  we  now 
turn  to  a  straightforward  linear  approach  for  modeling  the  effects  of  noise  and  predicting 
the  distributions  of  estimated  polarizabilities.  This  simple  method  makes  no  restrictive 
assumptions  about  the  distribution  of  the  noise  and  directly  reproduces  the  variability  of 
polarizabilities  observed  in  field  data. 
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4.1.2.  Predicting  the  distributions  of  polarizabilities.  Equation  23  relates  the  observed  noise 
e  to  the  perturbation  in  the  model  Am  via  the  forward  modeling  operator  G 

(23)  e  =  GAm. 

From  this  relation,  we  can  express  the  covariance  of  the  noise  as 

cov(e)  =  E(eeT) 

(24)  =  GF((Am)(Am)T)GT 

=  Gcov(m)GT. 

We  can  estimate  the  covariance  of  the  polarizabilities,  cov(m)  directly  from  the  set  of  esti¬ 
mated  target  polarizabilities  at  a  site,  as  displayed,  for  example,  in  figure  12.  Note  that  for 
each  target  in  this  set,  there  is  a  separate  G  corresponding  to  the  estimated  orientation  and 
location  of  that  target. 

We  now  assume  that  there  is  some  equivalent  dipole  source  at  the  median  location  req  of 
all  targets  in  our  observed  target  set.  For  example,  at  Pole  Mountain  and  Camp  Beale  the 
median  location  of  ISOs  is  directly  below  the  center  of  the  array,  at  approximately  20  cm 
depth.  The  forward  operator  Geq  at  this  location  can  then  be  used  to  compute  the  effective 
noise  covariance  using  equation  24. 

Geq  also  depends  on  the  orientation  of  the  dipole  source.  As  in  a  GSV  analysis,  we 
can  restrict  ourselves  to  targets  in  a  particular  orientation.  Alternatively,  we  can  specify 
a  probability  distribution  of  target  orientations  and  marginalize  over  azimuth  and  dip  to 
determine  the  expected  covariance  of  the  estimated  polarizabilities.  Appendix  B  develops 
expressions  for  the  covariance  of  the  estimated  polarizabilities  under  a  uniform  distribution 
of  target  dip  and  azimuth.  For  simplicity,  in  this  section  we  consider  only  vertically-oriented 
targets. 

Figure  20  compares  the  effective  noise  covariances  for  Pole  Mountain  and  Camp  Beale 
computed  using  this  approach.  Consistent  with  our  characterizations  of  the  noise  in  the 
previous  section,  there  is  significant  off-diagonal  structure  in  figure  20(a),  indicating  the 
presence  of  correlated  noise  on  the  data.  The  estimated  standard  deviations  in  figure  20(b) 
show  that  a  single  standard  deviation  at  each  time  channel,  as  assumed  in  the  previous 
section,  is  not  an  accurate  characterization  of  the  noise.  Instead,  we  find  that  are  three 
standard  deviations  required  at  each  channel:  one  corresponding  to  maximum  coupling 
between  transmitter  and  receiver  (e.g.  x  transmitter  and  x  receiver  component),  and  two 
components  for  perpendicular  transmitter/receiver  combinations. 

Inverting  at  a  fixed  target  location  is  a  linear  transformation  of  the  data,  and  so  with  the 
noise  covariance  from  equation  24,  we  can  directly  compute  the  covariance  of  the  polariz¬ 
abilities  at  the  location  rq  with  corresponding  forward  operator  G j  as 

cov(m)j  =  Gjcov(e)GjT, 


(25) 
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Figure  20.  Structure  of  the  noise  covariance  matrices  at  Camp  Beale  and 
Pole  Mountain,  (a)  Correlation  matrix  structure.  Each  block  in  the  correlation 
matrix  is  comprised  of  the  21  measurements  (7  receivers  x  3  transmitters) 
made  with  a  given  receiver  component  (b)  Estimated  standard  deviation  of 
the  noise  for  center  MetalMapper  receiver  (all  components)  at  Camp  Beale 
and  Pole  Mountain. 


with  G]  denoting  the  pseudo-inverse.  For  r.;  =  req ,  we  find  cov(m)?  =  cov(m)e9.  This  means 
that  a  target  at  our  equivalent  source  location  will  produce  a  distribution  of  polarizabilities 
that  is  exactly  that  which  we  recovered  from  the  field  data. 

Now  as  we  vary  the  target  location  r^,  there  is  a  commensurate  change  in  our  model 
uncertainty,  as  shown  in  figure  21.  This  analysis  is  still  restricted  to  vertically  oriented  tar¬ 
gets.  The  change  in  model  uncertainty  reflects  the  relative  change  in  curvature  of  the  misfit 
function  as  the  target  location  varies.  Figure  22  shows  the  dependence  of  the  polarizability 
uncertainty  on  the  horizontal  location  of  the  target.  The  uncertainty  is  relatively  constant 
within  the  sensor  footprint. 

At  each  possible  location  for  a  target,  we  now  have  the  machinery  to  compute  the  uncer¬ 
tainty  in  the  polarizabilities  that  is  expected  given  a  specified  noise  covariance.  In  practice, 
we  have  a  sample  of  targets  S  at  different  locations  and  wish  to  compute  the  overall  covari¬ 
ance  of  this  sample.  It  is  straightforward  to  show  that  the  covariance  of  the  sample  can  be 
estimated  as  the  mean  of  the  covariances  of  all  targets  in  the  sample 

M 


(26) 


cov(m)s  = 


M 


More  generally,  for  a  specified  distribution  of  target  locations  p( r)  the  model  covariance 
is  computed  via  the  expectation 


(27) 


cov(m)  =  y'cov(m)(r)p(r)(ir. 
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Figure  21.  Effect  of  target  depth  on  polarizability  uncertainty,  (a)  Predicted 
standard  deviation  of  polarizabilities  as  a  function  of  target  depth.  Predictions 
are  for  channel  one  MetalMapper  data  at  Pole  Mountain,  (b)  ISO  reference 
polarizability  (solid  lines)  and  predicted  95%  confidence  interval  (dashed  lines) 
for  ISOs  at  0.1  m  depth,  (c)  As  in  (b),  at  0.4  m  depth. 

For  specified  horizontal  and  vertical  spatial  distributions  of  targets,  we  approximate  the 
above  integral  at  a  discrete  set  of  points,  with  each  point  weighted  by 

(28)  P(r,)  =  =# 

Z^k  =  1 

Figure  23  shows  the  effect  of  assuming  normal  or  uniform  horizontal  distributions  of  ISO 
targets  on  the  predicted  polarizability  distributions.  The  uniform  horizontal  distribution  in 
(a)  might  represent  a  scenario  where  the  operator  of  the  sensor  platform  does  a  relatively 
poor  job  of  positioning  the  sensor  on  top  of  a  target.  If  we  replace  this  operator  with  a 
more  experienced/careful  driver  in  (b),  there  is  a  moderate  expected  reduction  in  polariz¬ 
ability  variance.  In  cases  (a)  and  (b),  we  have  assumed  a  uniform  vertical  distribution  of 
ISOs  down  to  a  maximum  depth  of  40  cm1.  Consequently,  the  overall,  averaged  variance 


1This  maximum  depth  may  correspond  to  the  clearance  depth  specified  by  a  regulator. 
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Figure  22.  Dependence  of  polarizability  standard  deviations  on  horizontal 
target  location,  Gridded  images  show  log10-transformed  standard  deviations. 
Markers  indicate  MetalMapper  receiver  locations.  Calculations  are  for  a  target 
at  20  cm  depth,  using  data  covariance  derived  from  Pole  Mountain. 

from  equation  27  is  dominated  by  targets  at  depth  and  the  benefits  of  improved  horizontal 
positioning  are  attenuated.  In  (c)  and  (d)  we  now  reduce  the  maximum  clearance  depth  for 
ISOs  to  20  cm.  This  has  a  strong  effect  on  the  polarizability  variance  for  both  drivers,  with 
a  dramatic  reduction  in  variance  for  our  experienced  driver  in  (d).  We  conclude  from  these 
simple  examples  that  the  important  controlling  variable  in  polarizability  variance  with  a 
cued  sensor  is  target  depth.  Good  horizontal  positioning  is  important,  but  may  only  provide 
modest  improvements  in  parameter  estimation  when  we  are  tasked  with  classifying  targets 
near  their  maximum  detection  depth. 

Thus  far  we  have  focused  exclusively  on  MetalMapper  data,  but  the  technique  is  read¬ 
ily  extended  to  any  cued  sensor,  with  the  functionality  to  compute  G  implemented  for  any 
data  type  we  already  invert  (TEMTADS2x2,  MPV,  etc.).  Figure  24  compares  recovered  ISO 
polarizability  distributions  for  the  BUDhh  and  TEMTADS2x2  sensors  deployed  at  Camp 
Beale.  The  BUDhh  measures  data  over  a  very  short  time  window,  with  the  result  that 
there  is  very  little  random  jitter  in  the  recovered  polarizabilities,  relative  to  the  late  time 
polarizabilities  measured  with  the  TEMTADS2x2.  Similar  to  its  larger  5x5  antecedent,  the 
noise  on  the  2x2  array  could  be  substantially  reduced  by  simply  averaging  the  polarizabil¬ 
ities  in  windows  of  e.g.  5  channels.  Evident  for  both  sensors  are  correlated  shifts  of  the 
polarizabilities,  in  this  case  likely  due  to  a  strong  background  response  at  the  site. 

We  use  these  polarizability  distributions  to  generate  noise  covariances  for  each  sensor,  and 
then  to  predict  the  distributions  of  polarizabilities  that  would  arise  for  each  sensor  for  ISOs 
distributed  uniformly  within  each  sensor’s  footprint.  We  again  consider  how  different  target 
distributions  with  depth  affect  polarizability  variance.  In  the  first  scenario  (middle  row  of 
figure  24),  we  represent  the  distribution  of  ISOs  at  Camp  Beale  as  a  normal  distribution 
with  mean  depth  15  cm  and  a  standard  deviation  of  5  cm.  This  roughly  corresponds  to  the 
actual  depth  distribution  of  ISOs  at  Beale.  Unsurprisingly,  this  model  closely  reproduces 
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Figure  23.  Predicted  ISO  polarizabilities  for  four  spatial  distribution  scenar¬ 
ios  with  Pole  Mountain  noise  model,  (a)  Uniform  horizontal  distribution  across 
the  sensor  footprint,  uniform  vertical  distribution  between  0.05  and  0.4  m.  (b) 
Normal  distribution  truncated  on  the  horizontal  interval  x,  y  G  [—0.3  0.3]  m, 
with  standard  deviation  of  0.1  m  in  both  x  and  y.  Uniform  vertical  distribu¬ 
tion  between  0.05  and  0.4  m  depth,  (c)  As  in  (a),  but  with  maximum  vertical 
depth  constrained  to  0.2  m.  (d)  As  in  (b),  but  with  maximum  vertical  depth 
constrained  to  0.2  m. 


the  observed  polarizability  distributions  for  both  sensors.  Now  if  we  make  the  parameter 
estimation  problem  a  little  bit  harder  by  distributing  targets  uniformly  with  depth  down  to 
30  cm  (bottom  row  of  figure  24),  the  BUDhh  shows  a  significant  increase  in  polarizability 
variance,  relative  to  the  2x2.  This  is  likely  because  the  BUDhh  measures  fewer  data  at  each 
sounding  (30  data  at  each  channel)  than  the  2x2  (48  data  at  each  channel).  This  produces 
a  larger  rate  of  increase  in  the  polarizability  variance  with  depth  for  the  BUDhh  than  the 
TEMTADS2x2. 

In  the  preceding  discussion  we  have  used  a  straightforward  linear  analysis  to  develop 
methods  for  predicting  the  distribution  of  estimated  polarizabilities  given  site  specific  noise. 
Beyond  its  simplicity,  the  approach  is  appealing  in  that  it  combines  the  problem  physics  - 


Polarizabilities  Polarizabilities  Polarizabilities 


32 


TEMTADS2x2 


BUDhh 


101 

10-’ 


101 

10° 

10-1 


10° 

10-’ 


Time  (ms) 


Time  (ms) 


Figure  24.  ISO  polarizabilities  for  TEMTADS2x2  and  BUDhh  handheld 
sensors  at  Camp  Beale.  Top  row:  actual  ISO  polarizabilities  for  each  sensor 
at  Camp  Beale.  Middle  row:  Predicted  ISO  polarizability  distributions  for 
each  sensor  with  a  uniform  horizontal  distribution  of  targets  within  the  sensor 
foorprint.  Targets  are  normally  distributed  in  depth,  with  a  mean  depth  of  15 
cm,  and  5  cm  standard  deviation.  Bottom  row:  As  in  middle  row,  but  with 
targets  uniformly  distributed  with  depth  down  to  30  cm. 


encapsulated  in  the  forward  modeling  matrix  G  -  with  a  realistic  noise  model  described  by 
a  dense  noise  covariance. 
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In  a  linear  uncertainty  analysis  we  typically  assume  perfect  recovery  of  target  location. 
There  is,  however,  some  accounting  for  imperfect  location  recovery,  and  other  nonlinear 
effects,  implicit  in  our  estimation  of  the  noise  covariance.  This  is  because  we  use  a  set  of 
actual,  estimated  polarizabilities  for  a  target  class  to  back  out  the  noise  covariance.  Cor¬ 
related  shifts  of  these  recovered  polarizabilities  resulting  from  errors  in  location  estimation 
are  therefore  propagated  through  the  analysis.  We  can  further  quantify  the  magnitude  of 
uncertainty  in  position  estimation  at  each  site  via  a  linearized  apparaisal.  This  is  identical 
to  the  linear  uncertainty  analysis,  except  the  forward  matrix  G  in  equation  25  is  replaced 
by  the  sensitivity  (Jacobian)  matrix  J  of  first  derivatives  with  respect  to  target  position. 
Figure  25  shows  this  analysis  for  the  MetalMapper  at  Camp  Beale  and  Pole  Mountain.  At 
each  target  location  we  compute  the  standard  deviation  of  our  location  estimate  as 

(29)  ar  =  yj  trace  (covfr'j,  rj)) 

with  co v(ri,Tj)  denoting  the  covariance  between  ith  and  jth  elements  of  the  position  vector 
r.  In  figure  25(a)  and  (c)  we  show  the  uncertainty  in  target  position  at  30  cm  depth  for  a 
single  object  scenario.  Within  the  sensor  footprint  the  median  standard  deviation  in  target 
location  is  approximately  30  cm  and  15  cm  at  Camp  Beale  and  Pole  Mountain,  respectively. 
In  the  linearized  analysis  we  can  also  examine  the  effect  of  multi-object  scenarios  on  location 
uncertainty.  In  (b)  and  (d)  we  introduce  a  second  target  near  the  edge  of  the  array  and 
compute  the  standard  deviation  in  position  as  a  function  of  the  location  of  the  first  target.  If 
the  sensor  is  properly  centered  over  one  of  the  two  targets,  then  the  uncertainty  is  essentially 
unchanged.  It  is  only  when  the  two  targets  are  both  near  an  edge  of  the  array  that  we  expect 
a  large  uncertainty  in  recovered  target  location.  Improved  infield  QC  procedures  should 
prevent  this  scenario,  so  that  we  can  reasonably  expect  a  relatively  small  error  in  recovered 
target  location  with  cued  sensors. 

Thus  far,  our  analysis  is  tied  to  specific  sites  for  which  we  have  a  representative  sample  of 
ISO  items  to  characterize  model  variability.  Ultimately  our  goal  is  to  predict  classification 
performance  for  an  arbitrary  site  characterized  by  particular  conditions  (topography,  target 
density,  etc.)  If,  for  the  moment,  we  regard  Pole  Mountain  and  Camp  Beale  as  end  members 
in  the  spectrum  of  classification  difficulty,  then  we  might  envision  a  system  where  a  user 
might  select  a  reference  site  most  similar  to  their  problem,  and  then  use  the  approach 
developed  here  to  predict  polarizabilities  and  classification  performance.  As  more  data  sets 
with  seeded  ISOs  become  available,  we  will  obtain  a  more  diverse  set  of  site  conditions,  with 
some  even  more  challenging  that  Camp  Beale  (e.g.  the  planned  demonstration  in  Waikoloa, 
HI).  In  ongoing  work  we  will  work  to  develop  metrics  (target  density,  magnitude  of  soil 
response,  sensor  positioning,  etc.)  necessary  to  characterize  each  site.  We  will  also  explore 
ways  to  average  existing  noise  models  (i.e.  covariances)  to  make  predictions  for  novel  site 
conditions. 


y  (in)  y  (m) 
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Figure  25.  Standard  deviation  of  positional  error  as  a  function  of  target 
position.  (a),(c)  Single-object  scenario  at  Camp  Beale  and  Pole  Mountain, 
respectively.  (b),(d)  Two-object  scenario  at  Camp  Beale  and  Pole  Mountain, 
respectively.  Marker  indicates  fixed  location  of  second  target  in  top  right 
corner  of  the  array.  Gridded  image  then  shows  the  uncertainty  in  the  location 
estimate  for  the  first  target. 
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4.1.3.  Predicting  the  ROC.  The  next  step  in  predicting  classification  performance  is  to  use 
the  distributions  of  polarizabilities  to  compute  the  distribution  of  the  decision  statistic  for 
each  target  class.  At  this  point  we  make  the  assumption  that  estimated  polarizabilities  are 
lognormally  distributed,  with  variances  predicted  using  the  above  analysis.  This  ensures  that 
the  polarizability  distribution  is  non-negative.  The  log-transformed  polarizabilities  are  then, 
conveniently,  normally-distributed,  allowing  us  to  analytically  compute  the  distribution  of 
the  decision  statistic.  For  classification  with  multistatic  sensor  data,  the  decision  statistic  0 
is  typically  a  misfit  of  log-transformed  estimated  polarizabilities  Lest  with  respect  to  some 
reference,  or  library,  polarizability  Lref 

N 


(30) 


</>  =  -  log (L7‘)? 

1=1 
N 


with  Xt  =  log (Lfst)  —  log(L'e^j.  The  summation  extends  over  some  subset  of  the  polariz¬ 
abilities.  For  high  SNR  targets  we  often  use  all  three  principal  polarizabilities,  whereas  for 
noisier  cases  we  might  only  use  the  primary  polarizability  to  compute  the  decision  statistic 
and  rank  targets.  Other  similarity  measures  may  also  be  used  as  a  decision  statistic  for 
ranking  targets.  For  example,  we  often  use  a  heuristic  function  of  the  form 


(z,f‘)o  -  uyr 
l/NT.nfr 

i=i  / 

with  7  ^  0.1.  The  following  analysis  can  readily  be  applied  to  nonlinear  functions  of  this 
form;  we  can  use  a  first  order  approximation  to  compute  the  variance  of  the  polarizabilities 
under  arbitrary  nonlinear  transformations. 

Continuing  with  the  decision  statistic  in  equation  30,  we  require  the  expected  value  of 
the  log  transformed  polarizabilities  in  a  class  to  be  the  logarithm  of  the  true,  reference 
polarizabilities  of  that  class 


(31) 
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=  E 

1=1 


(32) 


£(iog(Lr*»  =  log  urh 


If  the  lognormally-distributed  polarizability  has  predicted  variance  var(L),  then  the  variance 
of  the  corresponding  log-transformed  variable  is  given  by 

1  +  y/l  +  4var(L)  exp(— 2£'(log(L)))\ 


(33) 


var(log(L))  =  log 
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The  above  expressions  can  be  used  to  compute  the  moments  of  the  polarizability  misfit 
for  a  specified  target  class  with  respect  to  an  arbitrary  reference  polarizability.  The  sum 
of  N  standardized  normal  random  variables  is  distributed  as  a  ;y2  random  variable  with 
N  degrees  of  freedom.  In  this  case,  the  X,t  are  not  standardized  (i.e.  normalized  by  their 
standard  deviation).  We  therefore  use  a  variant  of  the  central  limit  theorem:  the  summation 
of  independent,  but  not  identically-distributed,  variables  in  equation  30  will  tend  to  a  normal 
distribution  (Billingsley,  1995)  with  mean  and  variance  given  by  equation  34. 

As  an  example,  we  consider  discrimination  between  ISO  targets  and  a  clutter  item. 
Ground  truth  photos  are  shown  in  figure  26.  We  have  selected  a  non-TOI  item  from  Camp 
Beale  that  has  a  very  similar  primary  polarizability  to  an  ISO.  In  figure  27  we  simulate  the 
distributions  of  polarizabilities  that  would  be  obtained  for  these  items  assuming  the  Pole 
Mountain  noise  covariance,  and  a  uniform  spatial  distribution  for  both  targets  to  a  maximum 
depth  of  40  cm.  We  first  compute  the  decision  statistic  using  only  the  primary  polarizability 
(top  row  of  figure  27).  This  produces  some  overlap  in  the  distributions  of  0,  such  that  the 
resulting  ROC  requires  us  to  dig  approximately  5%  of  clutter  in  order  to  ensure  that  99% 
of  all  TOI  are  found.  If  we  instead  use  all  polarizabilities  for  classification  (bottom  row  of 
figure  27),  then  the  ROC  indicates  perfect  classification.  This  is  because  there  is  a  significant 
difference  between  transverse  polarizabilities  for  these  targets  at  the  specified  noise  level. 

Figure  28  considers  the  same  classification  problem  as  in  figure  27,  but  uses  the  noise 
covariance  derived  from  Camp  Beale.  The  increased  noise  at  this  site  degrades  classification 
performance  using  the  primary  polarizability,  but  perfect  classification  is  still  expected  when 
all  three  polarizabilities  are  used  to  rank  targets. 
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(a)  ISO 


(b)  Clutter,  Camp  Beale  target  2530 


Figure  26.  Targets  used  for  ROC  prediction. 


Polarizabilities 


10  20  30  40 


Decision  statistic 

0.8 
0.6 
0.4 
0.2 
0 

0  5  10  15 


ROC 


Figure  27.  ROC  prediction  for  the  MetalMapper,  Pole  Mountain  noise.  In 
this  scenario  we  discriminate  between  ISOs  and  a  clutter  item  with  a  similar 
primary  polarizability.  Top  row:  Classification  using  primary  polarizabilities 
only.  Bottom  row:  Classification  using  all  polarizabilities. 
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Figure  28.  ROC  prediction  for  the  MetalMapper,  Camp  Beale  noise.  Plots 
are  the  same  as  in  figure  27,  only  the  noise  covariance  is  changed. 
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Thus  far,  we  have  restricted  our  performance  prediction  analysis  to  targets  in  a  single 
orientation,  with  the  results  in  this  section  generated  for  horizontally-oriented  targets.  Fig¬ 
ure  29  compares  ROC  prediction  for  the  same  TOI  and  non-TOI  items  in  horizontal  and 
vertical  orientations,  as  well  as  for  a  uniform  angular  distribution  of  target  orientations. 
Expressions  for  this  latter  case  are  developed  in  appendix  B.  Horizontal  targets  are  the  most 
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Figure  29.  Effect  of  target  orientation  on  decision  statistic  and  ROC.  Top 
row:  decision  statistic  and  ROC  for  classification  with  primary  polarizabilities. 
Bottom  row:  decision  statistic  and  ROC  for  classification  with  all  polarizabil¬ 
ities. 

difficult  scenario  for  discrimination  using  the  primary  polarizabilities.  In  this  orientation,  the 
primary  polarizability  is  less  coupled  to  excitation  by  the  horizontal  MetalMapper  coil.  This 
is  alleviated  somewhat  by  the  (horizontal)  fields  transmitted  by  the  vertical  x  and  y  coils, 
but  since  the  center  of  these  loops  is  vertically  displaced  from  the  receiver  plane,  data  from 
the  vertical  transmitter  coils  is  somewhat  lower  SNR  than  from  the  horizontal  coil.  This 
results  in  more  variability  in  the  estimated  primary  polarizability  for  horizontally-oriented 
targets  than  for  vertically-oriented  targets.  We  verify  this  effect  with  numerical  simula¬ 
tions  in  figure  30.  We  conclude  that  horizontal  targets  represent  the  worst  case  scenario  for 
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both  detection  and  classification  (when  we  rely  exclusively  on  primary  polarizabilities  for 
classification) . 

Unsurprisingly,  a  uniform  angular  distribution  of  targets  gives  a  classification  result  that 
is  intermediate  to  purely  horizontal  and  vertical  cases.  The  distributions  of  the  decision 
statistic,  using  only  the  primary  polarizability,  are  quite  similar  for  horizontal  and  uniform 
cases.  This  is  expected:  a  uniform  distribution  with  respect  to  target  dip  is  biased  towards 
horizontal  orientations2. 

Moving  now  to  discrimination  with  all  polarizabilities  (bottom  row  of  figure  29),  the  per¬ 
formance  of  vertical  and  horizontal  cases  is  reversed.  Secondary  and  tertiary  polarizabilities 
are  better  constrained  in  horizontal  orientations,  and  so  inclusion  of  these  features  allows 
for  perfect  discrimination  between  the  two  targets  considered  in  this  example.  Conversely, 
it  is  harder  to  excite  the  transverse  response  of  a  vertical  target,  and  using  transverse  polar¬ 
izabilities  for  classification  of  vertical  targets  therefore  degrades  performance.  The  uniform 
angular  distribution  is  again  intermediate  to  vertical  and  horizontal  cases. 


Figure  30.  Numerical  simulations  of  primary  polarizability  estimation  for 
vertical  and  horizontal  ISO  targets.  Targets  are  distributed  on  a  uniform  grid 
within  the  sensor  footprint  down  to  a  maximum  depth  of  40  cm.  Independent 
Gaussian  noise  is  added  at  each  channel  using  standard  deviations  estimated 
from  Pole  Mountain  background  measurements.  The  synthetic  data  for  each 
location  and  noise  realization  are  then  inverted  to  produce  the  polarizability 
estimates. 


In  the  preceding  analyses,  we  have  considered  a  simple  classification  problem:  we  have 
one  set  of  TOI  polarizabilities  and  one  set  of  non-TOI  polarizabilities.  We  then  predict  the 
ROC  that  would  be  generated  for  a  given  spatial  and  angular  distribution  of  each  target.  By 
predicting  classification  performance  for  specific  items,  this  analysis  can  give  a  site  manager 
a  concrete  sense  of  the  classification  capabilities  of  a  sensor  under  prescribed  conditions. 

2  At  a  horizontal  dip  angle,  an  azimuthal  angle  g  subtends  a  longer  arc  than  at  near  vertical  angles,  so  that 
a  uniform  distribution  on  the  surface  of  a  sphere  produces  more  samples  at  the  equator  than  at  the  pole. 
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A  more  general  and  realistic  classification  problem  is  to  predict  classification  between  a 
number  of  TOI  classes  and  a  prescribed  set  of  non-TOI  items.  This  is  a  straightforward 
generalization  of  the  approach  developed  here:  we  compute  the  distribution  of  the  deci¬ 
sion  statistic  for  each  item,  and  then  represent  the  distribution  of  the  decision  statistic 
within  each  class  (TOI  and  non-TOI)  as  a  mixture  of  the  distributions  of  each  item  in  that 
class.  The  component  distributions  in  the  sum  can  be  weighted  by  probabilities  reflect¬ 
ing  the  relative  frequencies  of  each  item.  Alternatively,  we  can  treat  the  non-TOI  class 
as  fundamentally  uncertain  and  define  an  underlying  non-TOI  distribution.  The  predicted 
polarizability  distribution  arising  for  specified  parameters  (spatial  and  angular  distributions 
and  noise  covariance)  can  then  be  convolved  with  the  underlying  distribution  of  non-TOI 
polarizabilities. 

The  choice  of  representative  non-TOI  items  is  an  important  consideration  when  predicting 
classification  performance.  We  have  maintained  a  database  of  targets  from  past  ESTOP 
demonstrations,  and  we  envision  a  decision  support  system  with  access  to  this  database. 
Non-TOI  could  be  selected  based  on  visual  similarity  to  TOI,  or  based  on  the  similarity  of 
polarizabilities  to  reference  polarizabilities. 
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4.2.  Performance  prediction  for  dynamic  sensors.  In  order  for  the  performance  pre¬ 
diction  tools  developed  in  the  previous  section  to  be  useful  for  site  managers  and  field 
geophysicists,  a  similar  analysis  must  be  available  for  dynamic  sensors.  We  must  be  able  to 
objectively  show  that  the  additional  time  and  expense  spent  on  collecting  and  processing 
cued  interrogation  data  will  provide  tangible  improvements  in  classification  performance, 
and  an  overall  reduction  in  cost.  Our  primary  focus  for  this  analysis  is  on  EM-61  data, 
since  this  sensor  is  familiar  to  practicing  geophysicists  and  provides  a  baseline  performance 
metric.  The  approach  can  be  readily  generalized  to  other  dynamic  sensors,  including  arrays 
of  EM-61s  and  newer  dynamic  platforms. 

As  with  cued  sensors,  we  can  estimate  a  site-specific  noise  covariance  using  a  linear  trans¬ 
formation  of  the  polarizabilities  for  ISO  items.  Figure  31  shows  total  polarizabilities  for 
ISOs  estimated  from  EM-61  cart  data  acquired  at  Camp  Beale.  These  features  exemplify  a 


Figure  31.  Total  polarizabilities  estimated  for  ISO  targets  at  Camp  Beale. 

notorious  difficulty  with  monostatic  data  processing:  poor  constraints  on  target  depth  pro¬ 
duce  a  large  spread  in  the  amplitude  of  recovered  polarizabilities.  This  is  further  illustrated 
in  the  size-decay  feature  space  for  the  Camp  Beale  handheld  area  in  figure  32.  The  poor 
constraint  on  estimated  depth  is  manifested  by  the  spread  of  the  size  parameter  within  each 
target  class.  The  decay  parameter,  however,  is  relatively  well  constrained.  In  particular, 
small  ISOs  have  a  fairly  tight  distribution  with  respect  to  decay.  This  is  also  evident  in  fig¬ 
ure  31:  the  primary  polarizabilities  have  a  reasonably  consistent  slope.  We  have  found  that 
the  simplest  and  most  reliable  approach  for  classification  with  EM-61  data  is  a  threshold  on 
the  decay  parameter. 

Using  the  approach  developed  in  the  previous  section,  we  can  predict  the  distribution  of 
estimated  polarizabilities  for  a  given  spatial  distribution  of  targets.  In  this  case  we  only 
specify  the  vertical  distribution  of  targets,  the  horizontal  distribution  of  targets  will  always 
be  assumed  uniform  across  the  dynamic  sensor’s  swath.  However,  a  new  complication  arises 
when  we  consider  dynamic  data:  the  dimension  of  the  forward  matrix  G  depends  upon  the 
along-track  and  cross-track  data  density.  This  means  that  we  cannot  directly  compute  the 
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Figure  32.  EM-61  size-decay  feature  space,  Camp  Beale  handheld  area. 


effect  of  varying  data  density  on  polarizability  variance,  since  the  dimension  of  the  estimated 
data  covariance  at  a  given  density  will  not  match  the  dimension  of  G  for  a  different  density. 

To  address  this  problem,  we  first  associate  the  observed  covariance  of  the  polarizabilities 
cov(m)  within  a  target  class  with  the  forward  modeling  Gi.  Matrix  Gi  is  evaluated  at 
a  median  target  depth,  with  a  specified  target  orientation  (horizontal  or  vertical),  and  a 
nominal  along-track  and  cross-track  data  spacing.  The  model  covariance  at  this  location, 
without  scaling  by  the  data  covariance,  is 

(36)  Si  =  (Gf  Gi)  1. 

We  denote  the  diagonal  elements  of  Si  -  the  unsealed  variances  of  the  polarizabilities  -  as 
v\.  For  G2  evaluated  for  a  different  data  spacing  (at  the  same  target  depth  and  orientation), 
we  can  similarly  generate  v2.  We  then  compute  the  ratio  of  standard  deviations  at  the  two 
data  densities 

(37)  Pi  = 

Finally,  we  form  the  diagonal  matrix 

(38)  R  =  diag(p). 

and  scale  the  model  covariance  of  the  target  class  by  R  to  obtain  the  expected  model 
covariance  at  the  data  density  corresponding  to  G2 

(39)  cov(m)2  =  Rcov(m)RT. 

This  transformation  preserves  the  correlation  structure  of  cov(m)  and  scales  the  model 
covariance  by  the  change  in  polarizability  uncertainty  dictated  by  the  change  in  data  density. 
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Figure  33  illustrates  this  effect  for  varying  cross-track  line  spacing:  as  the  line  spacing  is 
increased,  the  uncertainty  in  the  total  polarizability  increases  nonlinearly. 


Figure  33.  Dependence  of  EM61  total  polarizability  standard  deviations  at 
Camp  Beale  on  cross  track  spacing.  Standard  deviations  are  displayed  at  all 
four  EM61  channels,  with  the  magnitude  of  the  uncertainty  decreasing  with 
time. 


We  now  turn  to  computation  of  the  distribution  of  the  decay  parameter.  The  total  polar¬ 
izability  at  each  time  channel  has  variance 

(40)  var(Lt0ta/(4))  =  E  cov(Lj(bc),  Lj(tk )). 

Treating  the  total  polarizability  as  a  lognormally  distributed  random  variable,  the  decay 
parameter 

(41)  decay (tk,  tj)  =  Ltatal^  , 

is  then  the  ratio  of  two  lognormally  distributed  variables,  which  we  now  denote  as  X/Y . 
Taking  the  log  of  the  ratio  we  have 

(42)  Z  =  log  (0  =  log(X)  -  log(F). 

The  log  transformations  of  X  and  Y  are  normally  distributed,  so  that  Z  is  itself  normally 
distributed  with  mean  and  variance 

E(Z)  =  E(  log(X))  -  El  (log  (T)) 

(43) 

var(Z)  =  var(log(X))  +  var(log(T))  —  2cov(log(X),  log(F)). 

The  ratio  X/Y  =  exp(Z)  is  then  lognormally  distributed  with  mean  and  variance 

E(X/Y )  =  exp  (E(Z)  +  var  (Z)) 

(44) 

var  (X/Y)  =  (exp(var(E))  —  1)  exp(2  E(Z)  +  var  (Z)). 
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The  predicted  distribution  of  the  decay  parameter  depends  on  the  variance  of  the  total 
polarizabilities  at  the  selected  channels  tk  and  t:r  For  EM61  data  the  decay  parameter  is 
usually  computed  at  first  (t3  =  0.216  ms)  and  fourth  (tk  =  1.225  ms)  time  channels.  We  also 
find  that  the  correlation  between  the  total  polarizabilities,  quantified  by  the  covariance  term 
in  equation  43,  affects  the  decay  distribution.  Finally,  the  above  transformations  introduce 
a  dependence  on  the  size  of  the  target,  so  that  two  items  with  the  same  expected  decay 
parameter  but  different  total  polarizability  amplitudes  will  produce  different  distributions 
of  the  decay  parameter.  These  effects  are  illustrated  in  figure  34. 


(a)  Probability  density 


(b)  Standard  deviation 


(c)  Mean 


Figure  34.  Dependence  of  lognormal  decay  distribution  on  target  size  and 
correlation  coefficient.  Top  row  and  bottom  rows  are  for  correlations  of 
R  =  0.99  and  R  =  0.95,  respectively,  (a)  Dependence  of  lognormal  decay 
distribution  on  target  size.  As  target  size  decreases,  the  mode  of  the  decay 
parameter  distribution  shifts  to  smaller  values,  (b)  Standard  deviation  of  de¬ 
cay  parameter  as  a  function  of  size  scaling  parameter  a.  (b)  Expected  value 
of  the  decay  parameter  as  a  function  of  size  scaling  parameter  a.  Dashed  line 
indicates  the  true  value  of  the  decay  parameter. 


Taking  the  median  ISO  total  polarizability  at  Camp  Beale  as  our  reference,  in  figure  34 
we  scale  the  total  polarizability  by  a  factor  a  ranging  between  10-5  and  101,  leaving  the  true 
value  of  the  decay  parameter  (decay  ~  0.18)  unaffected.  In  (a)  we  show  the  dependence 
of  the  decay  distribution  on  this  scaling  of  target  size.  As  a  is  decreased,  the  distribution 
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becomes  increasingly  skewed,  with  the  mode  of  the  distribution  shifting  to  smaller  values.  In 
(b)  and  (c)  we  show  that  a  decrease  in  target  size  (smaller  a)  produces  an  increased  variance 
in  the  decay  parameter,  and  an  increased  positive  bias  in  the  expected  value.  These  effects 
are  consistent  with  what  we  observe  in  held  data.  Smaller  targets  produce  lower  SNR  data, 
and  so  we  expect  an  increase  in  variance  with  decreased  size.  We  also  find  that  at  later 
times  the  signal  falls  below  the  noise  floor  for  small  targets.  This  has  the  effect  of  biasing 
the  estimate  of  the  decay  parameter  upwards,  as  predicted  in  this  analysis. 

The  top  and  bottom  rows  in  figure  34  show  the  effect  of  varying  the  correlation  coefficient 
between  total  polarizabilities  at  channels  tj  and  tk-  Increasing  the  correlation  effectively 
decreases  the  noise,  and  reduces  the  variance  and  bias  of  the  decay  distribution.  For  the 
Camp  Beale  ISOs  the  correlation  coefficient  is  R  —  1,  so  that  the  expected  value  of  the 
decay  distribution  is  independent  of  target  size.  In  figure  34  we  reduce  the  correlation 
coefficient  to  R  =  0.99  and  R  =  0.95.  This  produces  a  positive  bias  in  the  expected  value. 
When  predicting  performance  for  new  data  sets,  the  correlation  coefficient  R  estimated  from 
known  ISOs  will  be  used. 

Figure  35  shows  the  predicted  distributions  of  the  EM61  decay  parameter  for  the  same 
targets  considered  previously  for  cued  sensor  performance  prediction  (see  figure  26).  Both 
targets  are  uniformly  distributed  down  to  40  cm  depth.  For  this  example,  the  clutter  item 
is  slightly  slower-decaying  than  an  ISO,  and  so  has  a  larger  expected  decay  parameter.  This 
is  the  opposite  of  the  usual  situation;  clutter  tends  to  be  faster  decaying).  However,  we  can 
still  predict  classification  performance  -  in  this  case  we  threshold  from  smallest  to  largest 
decay  parameter.  The  resulting,  mediocre  ROC  is  typical  of  EM-61  performance,  with  a 
false  alarm  rate  of  0.51  required  to  find  99%  of  the  ISOs. 
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Decision  statistic  ROC 


Figure  35.  Performance  prediction  for  the  EM61,  Camp  Beale  noise.  Left: 
distribution  of  the  decay  parameter  for  ISOs  and  a  slow-decaying  clutter  item 
(see  Camp  Beale  target  2530,  figure  26).  Both  targets  are  uniformly  dis¬ 
tributed  down  to  a  maximum  depth  of  40  cm.  Right:  predicted  receiver 
operating  characteristic. 
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5.  Conclusions  and  further  work 

In  this  report  we  initially  studied  the  prediction  of  target  response  thresholds.  Building 
on  previous  GSV  analyses,  we  showed  that  the  worst  case  detection  threshold  at  a  specified 
clearance  depth  depends  on  sensor  geometry  and  does  not  always  correspond  to  a  cross¬ 
track  azimuthal  orientation.  In  addition,  we  modeled  the  statistics  of  anomaly  amplitude 
for  a  uniform  distribution  of  cross-track  location  and  azimuthal  orientation.  We  found  that 
raising  the  detection  threshold  slightly  above  the  worst  case  scenario  can  drastically  reduce 
the  number  of  target  picks,  while  still  maintaining  a  high  probability  of  detecting  all  targets 
of  interest  at  the  maximum  clearance  depth.  Further  work  will  extend  this  analysis  to  newer 
detection  sensors  such  as  the  OPTEMA  and  dynamic  MetalMapper.  We  will  also  examine 
how  the  choice  of  time  channels  can  affect  the  number  of  target  picks  in  detection  data.  We 
are  preparing  a  journal  article  on  this  extended  GSV  analysis  for  publication  in  the  next 
year. 

Our  work  on  performance  prediction  focused  on  developing  realistic  models  of  the  noise  at 
each  site.  We  use  existing  data  sets  to  gain  an  understanding  of  how  recovered  polarizabilities 
can  vary  under  realistic  noise  conditions.  The  physical  uniformity  of  seeded  ISO  items  allows 
us  to  isolate  the  effects  of  site  specific  noise  on  the  polarizability  distributions.  We  then  use 
the  sample  covariance  of  the  ISO  polarizabilities  to  determine  the  corresponding  covariance 
of  the  noise  for  a  target  at  the  median  estimated  location  of  ISOs  in  the  sample.  This  noise 
covariance  is  highly  correlated  and  does  not  assume  a  constant  standard  deviation  for  all 
data  at  a  single  channel.  With  this  noise  covariance  we  can  then  predict  the  distribution  of 
polarizabilities  and  ROC  that  would  arise  for  any  specified  spatial  and  angular  distributions 
of  targets.  Relative  to  conventional  Monte  Carlo  simulation,  this  approach  to  performance 
prediction  is  very  fast.  Indeed,  once  the  forward  modeling  matrix  is  pre-computed  on  a  grid 
of  locations,  we  can  instantaneously  predict  performance  for  any  targets. 

To  this  point  we  have  not  tied  the  estimated  covariances  of  ISO  polarizabilities  at  ESTCP 
demonstration  sites  to  specific  environmental  conditions  at  those  sites.  Directly  linking  in¬ 
puts  such  as  soil  conditions  and  target  density  to  polarizability  variance  will  be  a  focus  of 
ongoing  work.  This  will  allow  us  to  predict  performance  under  arbitrary  conditions,  rather 
than  restricting  ourselves  to  sites  that  are,  for  example,  “Beale-like.”  Monte  Carlo  simula¬ 
tions  can  be  used  to  develop  these  relationships.  For  example,  we  recently  simulated  the 
distributions  of  polarizabilities  that  will  arise  in  the  presence  of  strongly  magnetic  Hawaiian 
soils.  These  results  will  be  useful  for  developing  an  efficient  predictive  model  quantifying 
the  connection  between  soil  response  and  classification  performance. 
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Appendix  A.  Summary  of  ESTCP  demonstration  data  sets 


Sensors 

Targets  of  Interest 

Results 

Demonstration  Site:  Camp  Sibert 

•  Geonics  EM-61  cart 

•  MTADS  EM- 61  array 

•  MTADS  mag  array 

•  Cued  EM-63 

•  4.2”  mortars 

The  results  for  all  sensor  combinations  were  excellent.  When  inverted  coopera¬ 
tively  with  magnetics  data,  the  EM63  cued  interrogation  was  the  most  effective 
discriminator.  All  33  UXO  were  recovered  prior  to  25  false  alarms. 

The  results  from  the  EM-61  cart  were  also  very  good,  although  24  false-positives 
were  required  to  excavate  all  105  UXO. 

Demonstration  Site:  San  Luis  Obispo 

•  MTADS  mag  and  EM- 
61  arrays 

•  Geonics  EM61 

•  MSEMS 

•  TEMTADS 

•  MetalMapper 

•  BUD 

•  60  mm  mortar 

•  81  mm  mortar 

•  4.2”  mortars 

•  2.36”  rockets 

•  one  each  of  37  mm,  3” 
and  5”  projectiles 

Magnetometer  detection  and  discrimination  performance  at  this  site  was  quite 
poor.  For  EM-61  production  data,  the  time-decay  rate  estimated  from  the 
recovered  polarizabilities  provided  an  effective  ranking  scheme. 

For  the  TEMTADS  data  the  library  method  was  the  most  effective  with  204  of 
206  TOI  recovered  along  with  131  of  1076  non-TOI.  The  other  discrimination 
methods  were  also  effective,  generating  between  2  to  4  false-negatives. 

The  library  method  was  again  most  effective  for  the  MetalMapper  data  with 
the  excavation  of  203  of  204  TOI  and  175  of  1205  non-TOI. 

Demonstration  Site:  Camp  Butner 

•  Geonics  EM61  cart 

•  TEMTADS 

•  MetalMapper 

•  105  mm  projectile 

•  37  mm  projectile 

•  M48  fuze 

This  site  exhibited  significant  magnetic  soil  response  and  relatively  high  target 
densities.  EM-61  performance  was  quite  poor,  with  approximately  70  %  of 
non-TOI  excavated  in  order  to  find  all  TOI. 

Excellent  classification  performance  was  achieved  with  TEMTADS  data  pro¬ 
cessing,  all  TOI  were  readily  identified  and  only  4  %  of  clutter  were  dug.  The 
classification  method  relied  on  a  misfit  with  respect  to  all  polarizabilities,  fol¬ 
lowed  by  classification  with  total  polarizability  only. 

MetalMapper  performance  was  comparable  to  the  TEMTADS  for  identification 
of  90  %  of  TOI,  but  higher  noise  levels  in  one  of  the  sensors  deployed  at  this 
site  produced  2-3  missed  TOI  most  performers’  diglists. 

Sensors 

Targets  of  Interest 

Results 

Demonstration  Site:  Pole  Mountain 

•  Geonics  EM61  cart 

•  MetalMapper 

•  Stokes  mortar 

•  75  mm 

•  60  mm  mortar 

•  57  mm 

•  37  mm  projectile 

•  Small  ISO 

This  demonstration  was  divided  into  two  parts  (Years  1  and  2).  The  EM-61 
achieved  mediocre  performance  for  both  parts:  the  false  alarm  rates  (FAR) 
for  Year  1  and  Year  2  were,  64.0  %  and  68.0  %  respectively.  All  TOI  were 
recovered  at  the  specified  EM-61  operating  points. 

The  MetalMapper  data  produced  consistently  excellent  classification  perfor¬ 
mance.  For  example,  in  the  Year  1  test  the  Library  method  resulted  in  the 
excavation  of  all  TOI  with  a  FAR  of  2.77  %.  A  number  of  analysts  processed 
these  data  and  achieved  similar,  “near-perfect”  performance. 

Demonstration  Site:  Camp  Beale 

•  Geonics  EM61  cart 

•  MetalMapper 

•  MPV 

•  TEMTADS2x2 

•  BUDHH 

•  105  mm 

•  81  mm  mortar 

•  60  mm  mortar 

•  37  mm  projectile 

•  Small  ISO 

This  demonstration  tested  classification  in  a  portion  of  the  site  amenable  to 
vehicular  towed  systems  (i.e.  the  “Open  Area”).  In  addition,  man-portable 
sensor  data  were  collected  with  the  TEMTADS  2x2,  BUD,  and  MPV  sensors 
in  a  treed  section  of  the  site  (i.e.  the  “Portable  Area”). 

For  the  Portable  Area,  the  EM-61  analysis  required  approximately  30  %  of 
clutter  to  be  excavated  in  order  to  find  all  TOI.  For  the  Open  Area,  the  false 
alarm  rate  was  50  %  with  this  sensor. 

MetalMapper  data  were  acquired  by  Parsons  and  CH2M  HILL  in  the  Open 
Area,  with  a  different  analyst  processing  each  data  set.  Application  of  a  li¬ 
brary  classification  method  to  the  CH2M  HILL  data  resulted  in  99.2  %  of  TOI 
identified  (due  to  one  missed  TOI)  and  27.4  %  scrap  dug  at  the  stop  dig  point. 
For  the  Beale  Parsons  data  two  stage  (all  polarizabilities/total  polarizability) 
classification  with  a  support  vector  machine  was  able  to  identify  all  TOI  at  the 
stop-dig  point,  with  a  22  %  false  alarm  rate. 

All  portable  sensors  data  sets  produced  diglists  with  no  TOI  left  in  the  ground 
in  the  portable  area  of  the  site.  False  alarm  rates  were  also  quite  similar  for  all 
portable  sensors,  averaging  approximately  25  %  of  clutter  dug  in  order  to  find 

all  TOI. 
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Appendix  B.  Expected  polarizability  covariance  under  a  uniform 

DISTRIBUTION  OF  TARGET  ORIENTATIONS 


To  predict  classification  performance  in  section  4  we  use  a  single  target  orientation  (hor¬ 
izontal  or  vertical)  to  propagate  uncertainties  and  predict  the  distribution  of  the  decision 
statistic.  More  realistically,  we  can  model  a  specified  distribution  of  target  orientations.  A 
uniform  distribution  of  target  azimuth  and  dip  is  the  obvious  choice,  and  here  we  develop 
analytic  expressions  for  the  covariance  of  the  estimated  polarizabilities  under  a  uniform 
rotation. 

The  symmetric  polarizability  tensor  M  is  related  to  the  diagonal  matrix  of  principal 
polarizabilities  L  via  the  Euler  rotation  matrix  A 

(45)  M  (t)  =  AL(t)AT. 


For  target  dip  b  E  [0  7r]  and  azimuth  g  G  [0  27t],  we  define  A  as 


(46) 


A  = 


cos(g)  cos (b) 
—  sin (g)  cos (b) 
sin(6) 


sin(flr) 
cos (g) 
0 


—  cos(g)  sin(6)\ 
sin(g)  sin(fe) 
cos (b)  J 


Under  this  convention  the  diagonal  polarizability  matrix  L  is  ordered 


( L2(t ) 

0 

°  \ 

(47) 

m  = 

0 

Ls(t) 

0 

l  0 

0 

so  that  a  vertical  target  with  g  =  0  has  its  primary  polarizability  {L\ (t))  aligned  with  the  z 
axis.  The  predicted  data  d  for  a  target  at  location  r  are  given  by 


(48) 


d  =  G(r)m 


with  the  model  vector  m  comprised  of  the  unique  elements  of  the  polarizability  tensor 
(49) 


m 


/Mn\ 

M12 

Mis 

M22 

M23 

\M33/ 


/ 


L2  cos (b)2  cos(g)2  +  Li  cos (g)2  sin(6)2  +  L3  sin(5f)2  \ 

L2  cos (g)  sin(^)  cos(6)2  —  Li  cos (g)  sin (g)  sin(6)2  +  L3  cos(^)  sin(p) 
L2  cos(6)  cos(g)  sin (6)  —  L\  cos (6)  cos(g)  sin(6) 

L2  cos (b)2  sin(5f)2  +  L3  cos (g)2  +  L\  sin(5)2  sin(5f)2 
L\  cos (b)  sin(6)  sin (g)  —  L2  cos (6)  sin(6)  sin (g) 

Li  cos (b)2  +  L2  sin(6)2  / 


with  the  dependence  on  time  t  now  implicit.  In  this  formulation  the  forward  modeling  matrix 
G  has  dimensions  N  x  6,  with  N  the  number  of  data  measured  at  a  single  time  channel. 
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From  the  above  expressions,  we  define  the  operator  a 

(  cos  (6) 2  cos  (g) 2  sin(g)2 

—  cos(6)2  cos(g)  sin(g)  cos(g)  sin (g) 

cos (6)  cos (g)  sin(6)  0 

cos(6)2  sin(g)2  cos(g)2 

—  cos(6) sin(fi) sin(g)  0 

\  sin(6)2  0 


(50) 


a  = 


cos(g)2  sin(6)2 

—  cos (g)  sin(6)2  sin(g) 

—  cos(6)  cos (g)  sin(6) 

sm(b)2  sin(g)2 
cos (b)  sin(6)  sin(g) 
cos (b)2  ) 


such  that  a  maps  from  the  principal  polarizabilities  to  unique  elements  of  the  polarizability 
tensor 


(51) 


m  =  a 


(l2\ 

L3 

w 


Now  we  assume  that  the  covariance  of  the  data  for  a  given  target  orientation  is  given  by  a 
transformation  of  the  polarizability  covariance 


(52)  cov(d)  =  Geq  a  cov(L)  aT  G^q 

with  Geq  corresponding  to  the  median  location  (relative  to  the  sensor)  of  all  targets  used  to 
compute  cov(L).  Now  for  a  target  at  location  r  the  covariance  of  the  estimated  polarizabil¬ 
ities  can  be  expressed  as 

cov(L)  =/3  G^(r)  cov(d)  (G^(r))T  /3r 

(53) 

=f3  Gt(r)  Geq  a  cov(L)  aT  G^  (G^r))71  f3T 

where  the  pseudo-inverse  is  G^  =  (GTG)-1G.  The  operator  /3 

(  cos(6)2  cos(g)2  sin(g)2  cos(g)2  sin(6)2  \ 

— 2  cos(6)2  cos(g)  sin(g)  2cos(g)  sin(g)  —  2  cos(g)  sin(6)2  sin(g) 
cos(6)2  sin(g)2  0  — 2cos(6)  cos(p)  sin(6) 

2  cos (b)  cos (g)  sin(6)  cos(g)2  sin(6)2  sin(g)2 

2  cos(6)  sin(6)  sin(g)  0  2  cos (6)  sin (6)  sin^) 

sin(6)2  0  cos  (6)2 

corresponds  to  diagonalization  of  the  estimated  polarizability  tensor 

(55)  L  =  AtMA, 
so  that 

(56) 


(54) 


P  = 


V 


(l2\ 

u 

\U) 


=  /3m. 
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In  addition,  since  A  is  orthogonal,  we  find  that 

(57)  /3a  =  I. 

Equation  53  provides  an  expression  for  the  covariance  of  the  estimated  polarizabilities  at 
a  given  orientation.  We  now  assume  a  uniform  probablility  distribution  of  target  dip  b  and 
azimuth  g,  with  the  independent  densities  (Brannon,  2002) 

1 

27 r 

p(b)  =  -  sin  (b). 

We  marginalize  equation  53  over  azimuth  and  dip  to  obtain  the  expected  covariance  of  the 
polarizbility  estimates 

27 T  7 T 

(59)  E(cov(L))  =  —  J  J  cov(L)  sin (5)  db  dg 

g= 0  6=0 

These  integrals  involve  only  trigonometric  functions  and  can  be  evaluated  analytically.  The 
resulting  expressions  are  quite  lengthy  and  so  we  do  not  reproduce  them  here. 


